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ABSTRACT 

We combine a large set of quasar luminosity function (QLF) measurements from the rest-frame optical, soft 
and hard X-ray, and near- and mid-infrared bands to determine the bolometric QLF in the redshift interval 
Z = 0-6. Accounting for the observed distributions of quasar column densities and variation of spectral energy 
distribution (SED) shapes, and their dependence on luminosity, makes it possible to integrate the observations 
in a reliable manner and provides a baseline in redshift and luminosity larger than that of any individual survey. 
We infer the QLF break luminosity and faint-end slope out to z ~ 4.5 and confirm at high significance (> IOct) 
previous claims of a flattening in both the faint- and bright-end slopes with redshift. With the best-fit estimates 
of the column density distribution and quasar SED, which both depend on luminosity, a single bolometric 
QLF self-consistently reproduces the observed QLFs in all bands and at all redshifts for which we compile 
measurements. Ignoring this luminosity dependence does not yield a self-consistent bolometric QLF and there 
is no evidence for any additional dependence on redshift. We calculate the expected relic black hole mass 
function and mass density, cosmic X-ray background, and ionization rate as a function of redshift and find they 
are consistent with existing measurements. The peak in the total quasar luminosity density is well-constrained 
at z = 2.15 ±0.05. We provide a number of fitting functions to the bolometric QLF and its manifestations in 
various bands, and a script 3 to return the QLF at arbitrary frequency and redshift from these fits. 
Subject headings: quasars: general — galaxies: active — galaxies: evolution — galaxies: luminosity func- 
tion — cosmology: observations — X-rays: galaxies — infrared: galaxies — ultraviolet: 



galaxies 



1. INTRODUCTION 



Determining the nature and cosmological evolution of the 
quasar luminosity function (QLF) has been of inter est since 
quasa rs were first identified as cosmological sources ISchmidtl 
1968), and understanding the QLF is crucial to inferring the 
formation history of supermassive black holes, as well as the 
buildup of cosmic X-ray and infrared (IR) backgrounds and 
the contribution of quasars to reionization. Furthermore, the 
recognition that blac k holes appear to reside at th e centers 
of most galaxies (e.g., Kormendv & Richstone 1995) and that 
the masse s of these black holes are correlated with either 
the ma ss (Magorrian et alJ 119981) or the vel ocity dispersion 
jFerrarese & Merrittl2000tlGebhardt et all20 00) of their host 
spheroids, demonstrates a link between the origin of galaxies 
and supermassive black holes. Hence, determining the evo- 
lution of the QLF is also critical to understanding galaxy for- 
mation and evolution. 
The study of the QLF has a long history (e.g. , 



Schmidt & GreerJl983tlKoo & Kronll98llBoyle et alJl98£ 
Hewett et alJ Il993t iHartwick & Schadd 119901: IWarren et alJ 
19941 ISchmidt et alJl 19951 iKennefick et alii 19951 iPeill 19951) . 

but in recent years, surveys such as the Two Degree Field 
(2dF) QSO Redshift Survey (2QZ; iBovle et alJ 1200(3) and 
the Sloan Digital Sky Survey (SDSS: lYork"et alJl2000l) have 
provided large, homogeneous quasar samples ove r the range 
of redshifts z = - 6 (e.g., Boyle et al. 2000; Fan et al. 



2001a, 2004; Croom et al. 2004; Richards et al. 2005, 2006b; 


Jiang et alJ 




. In addition, a great deal of information 



on the X-ray and infrared properties of quasars has become 
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available, and surveys with e.g. Chandra, XMM, ROSAT, and 
Spitzer have enabled studies of the evolution of the QLF 
across many frequ encies (e.g..lMiv aii et alJ l2000t lUeda et all 
| 2003[ lHaas et all 1 20041: IBarger et alJ 120051 lHasinger et alJ 
l2005UBrown et alJ2006tlMatute et alJ2 006). 

Surprising and suggestive trends have emerged 
from these stu dies. For e xampl e, both the s pectral 
shapes (s-g - IWilkese talJ Il9 94t iGre enetal] 1199 
Vign aliet al] 120031 [Stratevaet alT 120051: iRichardset " 
2006c; Steff en et alJ 120061) and column densif 
distributions (e.g., iHill. Goodrich. & DePovl [199 

119991 IWillott et alJ l200*i 



Simpson. Rawlings. & Lacv 



Simpson & Rawlings 2000; Steffen et al. 2003; Uedaetal. 
2003; Grimes, Rawlings, & Willott 2004; Hasinge ri l2t„- ,, 
Sazonov & Revnivtseyl l2004t IBarger et alJ 120051: ISimpsoiJ 
20051 lHao et al] 120051) of quasars appear to depend system- 
atically on luminosity, with the brightest quasars being the 
least obscured and the most dominated by the optical/UV 
portion of the spectrum. X-ray studies (e.g., Page et al. 1997; 
Mivaii et alJl2000t _l200lt lLa Franca et al J 120021: ICowie et all 
20031 lUeda et alJ 120031; iFioreetalJ l2003t IBarger etal J 
2003b; Hasinger et al. 2005) an d som e optical, radio, and IR 
measurement s (e.g .. iHunt etalll2004t ICirasuolo et alJrlOOSl 
MatuteetaD[2006) indicate that the space density of low 
luminosity active galactic nuclei (AGN) peaks at redshifts 
lower than that of bright quasars, specifically observing a 
flattening of the QLF faint-end slope with redshift. It has 
been argued that this follows a similar pattern of "cosmic 
downsizing" as ha s recently been ob served in galaxy spheroid 
populations (e.g.. ICowie et alJll996l) . Similarly, the bright- 
end slope of the QLF appears to become shallowe r towards 
higher redshift s, from both d irect measurements iFan et alJ 
I2001bl 120031: iRichards et all l2006bl) and (albeit w eaker) 
constraints from gravi tational lensing ( Comerford et alj 2002t 
iWvithe & Loebl2002l:lWvithel2004tlRic"hards et alJ2006at) . 
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Various models have been proposed to explain these trends, 
many of which postulate that feedback from black hole 
growth plays a key role in determinin g the black hole- 
host galaxy (e.g. Mrh - a) relationships ( Silk & Reesl ll998t 
iDi Matteo et alJl2005l) . and co-evolution of black holes and 
their host spheroids. The evolution of the faint and bright- 
end slopes may be linked to these processes, with AGN feed- 
back providing the mechanism for cosmic downsizing, shut- 
ting down the growth of the most massive syste ms at high red- 
shift (e g..lMerlonill2004t iGranato et"aT]l2004t iHopkins etafl 
2006b; Croton et al. 2006) and potentially steepening the low- 
redsh ift bright-end QLF slope as a result ( Scannapieco & Oh 
l200l . while this feedback-driven quas ar decay d etermines 
the shape of the faint-end QLF (Hopkins et al. 2006a). 
Feedback may also explain trends in obscuratio n with lu- 
mino sity, either through dust sublimation (e.g., Lawrence 
1 19911) or expulsion o f gas and dust on galactic scales (e.g., 
Hopki ns et al.l l2005dh . But attempts to quantitatively link 
the downsi zing of quasar and galaxy populations, both the- 
oretically (Hopkins et alJ 12006 b) and observationally (e.g., 
IHopkins et al.ll2006dl) . depend on a reliable determination of 
the QLF. 

However, inferences drawn from the observed trends suf- 
fer from complications arising from various biases. For 
example, optical quasar surveys generally probe large vol- 
umes (~ 10,000deg 2 ), enabling uniform sample selection at 
many redshifts and the discovery of the brightes t quasars, but 
miss s ubstantial populati o ns of obscured (e.g., iRisaliti et alJ 
119991 lUeda et aj.l L2003t iTreister et al.l 120041) or heavily 
redde ned (e.g. Iwebster et alJ 119951 [Richards et all 1200 It 
Broth erton et all 120011 iGreggetaLl 120021 IHopkins et all 
2004) quasars, and generally cover a relatively small base- 
line (~ 1 -2dex) in luminosity. Several seminal efforts (upon 
which we seek to expand) have been made to extend these 
basel ines by compiling various optical quasar ob servations 
(e.g.. |Hartwick & Schaddl990tlWarren et alll994l) . most no- 
tably |Pej| ([T995), but these have still been severely limited by 
obscuration/reddening and span only ~ 2-2.5 dex in luminos- 
ity, and furthermore the area and depth of quasar surveys have 
since increased by orders of magnitude. Hard X-ray and IR 
samples provide a more complete census of obscured quasars, 
and sample much fainter luminosities and wider (~ 4-5 dex) 
luminosity ranges, but are measured from much smaller sur- 
vey areas (< 1 deg 2 ). Without carefully accounting for the dif- 
ferential effects of obscuration, different spectral shapes, and 
selection effects across bands, it is not clear if trends observed 
primarily at a single frequency are physically meaningful or 
robust across frequencies. Furthermore, theoretical models 
typically do not deal directly with the quasar luminosity in a 
given band, but instead treat the bolometric quasar luminos- 
ity, and it is not clear that a simple bolometric correction can 
reliably translate between the two — the effects above may 
change the observed QLF shape as a function of frequency, 
luminosity, and redshift. 

In this paper, we combine recent measurements of the 
QLF in many wavelengths from the mid-IR through hard 
X-rays, to determine the observed bolometric quasar lu- 
minosity function. By utilizing multiwaveleng th measure- 
ments of quasar SE Ds (lElvis et aljfl994l IVanden Berk et al. 
200!| iTelfer et alJ 12002b IVignali et alJ l2003t iFanetalJ 
20031 iHat ziminao g lou et all l2005t iRichards et alJ 12006c ; 
Shemmeretal. 2006) that probe the distribution of obscu- 
rati on up to and including Compton-t hick column densi- 
ties JRisaliti et alJl999Hl Jeda et ali2003UTreister et aH200i 



iMainieri et afll200l lHao^tai]|200l Fibzzi et alJl2006h . we 
can predict the QLFs that would be observed as a function 
of luminosity, frequency, and redshift from a given bolomet- 
ric QLF and compare these simultaneously with all observed 
QLFs. This makes it possible to constrain the shape and evo- 
lution of the QLF over ~ 8 orders of magnitude in luminosity 
and ~ 9 in space density, from z = 0-6, larger than the base- 
lines probed by any individual survey. 

In § |3 we describe the observational data sets , including 
the observe d co lumn density distributions (§ 12. 21 and SEDs 
adopted (§ 12. It . and consider in detail the consistency of 
the QLF across various frequenc ies g iven different simpli- 
fications of these distributions (§ 12.31 . In §[5] we calculate 
the bolometric QLF as a function of luminosity and redshift, 
and consider the detailed evolution of the QLF shape and its 
manifestation in different observed bands. In §0we con- 
sider a number of complementary constraints, testing models 
of quasar lifetimes, the evolution in t he QL F shape, and the 
buildup of the black hole population (§ 14. Il l, and we calculate 
th e reli c black hole mass functi on, c osmic X-ray background 
(§ I4.2> . and ionization rates (§ 14.3b expected from the bolo- 
metric QLF and compare these to observations. In § [5] we 
summarize and discuss our results and future prospects for 
improving the measurements. 

We adopt a n M = 0.3, fi A = 0.7, Hp = 70km s"' Mpc -1 cos- 
mology, consistent with Sper gel et alJ 12003. 2006). Quasar 
B-magnitudes are Vega, and L Q refers to the bolometric solar 
luminosity L Q = 3.9 x 10 33 ergs _1 . 

2. THE OBSERVATIONAL DATA SET 

In what follows, we compile a large number of binned QLF 
measurements in different redshifts ranges, observed wave- 
bands, and luminosity intervals. Table \l\ lists the samples 
used, with the survey and fields for each, the rest wavelength 
or band measured, the redshift and luminosity range of the 
observed quasars, and the number of quasars in the sample. 
We use the term "quasar" rather loosely throughout, as the 
traditional Mb < -23 cut is not readily applicable to multi- 
wavelength observations or obscured sources. Our compi- 
lation instead attempts to represent all AGN with intrinsic 
(obscuration-corrected) luminosities above the observational 
limits at each redshift. This generally extends to the typi- 
cally adopted ~ 10 42 ergs _1 X-ray luminosity limit (Lbd ~ 
10 10 Lq), below which confusion with normal star-forming 
and starburst galaxies becomes problematic (although re- 
solved hosts at low redshift extend this to ~ 1O 9 L0; e.g., 
lHao et alJ2005T) . 

Some of the observational estimates are derived from the 
same sets of observations, or different updates of the same 
data sets. To avoid "double-counting," at the few points where 
two binned QLFs overlap in both luminosity and redshift and 
are derived from the same sample (at that luminosity and 
redshift), we discard the less well-constrained (usually less 
recent) point. This is generally rare, and has a negligible 
effect on our results. It d oes reduce significantly the sam- 
ples of Mivaii et al. (20011) , for w hich the Type 1 QLFs ar e 
updated in lHasinger et al J ( 120051) . and ICroom et all ( 12004 . 
which h as many (espec i ally f aint) luminosity intervals up- 
dated in Richards et al. (2005). Some measured Q LFs do 
not appear in this compilation (e.g.. lBovle et alJ2"b 00). as up- 
dated versions from the same or expanded samples exist (e.g., 
ICroom et alJ2004h . We have re -fit our results with these sam- 
ples (as well as those of e.g. lPeilll995t Hia ng et alj[2006b ; 
ISazonov et aT1l2006t Isiiinozaki et al J 120061: iBeckmann et alJ 
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2006b) and find the differences are negligible. 

2.1. The SED and Bolometric Corrections 

We construct a model "intrinsic" (un-reddened) quasar 
SED to compare observations in different ba nds. The tem- 
plate s pectrum initially follows that derived in M arconi et alJ 
(2004) and consists of a broken power law in the optical- 
UV with q n = -0 44 (L„ oc u a °) for 1/xm < A < 1300 A 
( Vand en Berk et alJ 12001ft , and a VY = -1.76 from 1200- 
500 A (Telferetal. 2002), essentially the modal spectral 
slopes for optically bright blue quasars ( Richards et al. 2006c, 
appropriate given that this is supposed to be a pre-reddened 
spectrum). Because a number of slightly different bands have 
been used in optical quasar surveys (e.g. B, g, i, 1450 A; 
see Ta ble [TJ, we furth er overlay the template spectrum of 
iRichards et alJ ll2006d) derived from the optically blue (i.e. 
un-reddened) subsample onto this broken power-law. Techni- 
cally, we factor out the best-fit power-law in the given range 
from the Richar ds' et al.l (l2006d) spectrum and then multiply 
our template by the residuals, but t his mak es little difference 
compared to adopting the IRichards et alJ J2006d) spectrum 
over thi s range. We do th e same with the template spec- 
trum of Vande n Berk et al.l (|2001) for optical wavelengths 
1 fxm < A < 1300 A, if a detailed correction from an observed 
band to a more "standard" band is needed, but such correc- 
tions have generally been provided for the QLF measurements 
in Tabled 

Longwards of A > 1 um, we adopt the mean spectrum from 
Richa rds' et alJ (l2006d) . with a typical observed IR "bump" 
from reprocessing, eventually truncated as a Rayleigh-Jeans 
tail of blackbody emission (a = 2), and a no-obscuration zero 
point taken from the optically blue (un-reddened) subsam- 
ple template spectrum. This gives a 15/im to /?-band correc- 
tion similar to the typically ado pted log (Li<;/Lp) = 0.23 (e.g., 
Elvis et al]l 19941 iHat ziminao glou et alJl2005l Richards et alJ 
2006c:" lMatute et alJl2006t iBrown et alJ 120061) for optically 
unobscured (Type 1) quasars. 

The X-ray spectrum beyond 0.5 keV is determined by a 
power-law, with an intrinsic photon i n dex V = 1.8 (e.g. , 
Georg eetail H9981 iPerola et alJ 12001 iTozzi et all 12003) 
and an exponential cutoff at 500 keV. A refle ction com- 
ponent is included following JUeda et alJ (| 2003 | ). generated 
with the PEXRAV model (Magdziarz & Zdzi arskl 11995ft 
in the XSPEC package with a reflection solid angle of 
2tt, inclination cos(/) = 0.5 and solar abundances. The 
X-ray spectrum is then renormalized to a given a ox = 
-0.3841og[L y (2500A)/L !/ (2keV)], and the points at 500 A 
and 50 A are connecte d with a power-law . The value of a™ 
depends on luminosity (tWilkes et a ll 19941 iGreen et alll995l 
Vigna lTet al.l 12003*1 FStrateva et al. 2005), an d w e adopt the 
most recent determination by Steffen et al. (200a), 

a ox = -0.107 log(L,, 250 o/[ergs- I Hz- I ])+ 1.739, (1) 

determined specifically for unobscured (Type 1) quasars. The 
above equation is derived from the least-squares bisector of 
the Lj,(2500 A) — Lj,(2keV) relation, since neither luminosity 
can be properly taken as an "independent" variable - the dif- 
ference if e.g. L„(25 00A) is consid ered independent is gen- 
erally small (see e.g. ISteffen et al.ll2006T) . but can be of im- 
portance for the most luminous X-ray AGN (yielding a dif- 
ference of ~ 30% in the bolometric correction at Lo.5-2keV ~ 
10 47 ergs _I ). The baseline for these observations is suffi- 
ciently large that this relation has been determined for nearly 




9 10 11 12 13 14 
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FIG. 1 . — Bolometric corr ections for B-band, mid-IR, soft and hard X-ray 
bands, determined in ij |2.1l from a number of observations as a function of 
luminosity and given by the fitting formulae in Equation |2j- The lognormal 
dispersion in the distribution of bolometric corrections at fixed L, given by 
Equation |5J is shown as the shaded range for each band. The full quasar 
template spectrum is available for public download 3 . 



all luminosities of interest in the compiled QLF measure- 
ments. Recent comparisons between large sample s of quasars 
select ed by both optical and X-ray surveys ( iRisaliti & Elvisl 
120051) further suggests that this is an intrinsic correlation, 
not driven by e.g. the dependence of obscuration on lumi- 
nosity, as does our comparison of bolometric QLFs derived 
below. There is no evidence for a trend of a ox with red- 
shift, or for any other trend in spectral shape with re dshift 
(e.g. lElvis etallll99l IVanden Berk et^l200"il ITelfer et all 
1 20021: 1 Vignah et alJ2003HFan et all2003t ISteffen et alJ2006l: 
Richards et al. 2006c; Shemmer et al. 20061 but see also 
Bechtold et al. 2003), so our spectrum depends only on lu- 
minosity. 

Ultimately, the bolometric corrections (with zero attenu- 
ation) derived can be accurately approximated as a double 
power-law 




with (c u k u c 2 ,k 2 ) given by (6.25,-0.37,9.00,-0.012) 
for L band = L B , (7.40,-0.37,10.66,-0.014) for 
Li 5Aim , (17.87,0.28,10.03,-0.020) for L . 5 -2keV, and 
(10.83,0.28,6.08,-0.020) for L 2 -iokeV. The ki term is 
important when the given portion of the spectrum is not 
dominant, controlled by the scaling of a ox , and the k 2 « 
term represents the nearly constant bolometric correction 
when a given portion of the spectrum dominates the bolo- 
metric luminosity. Figure ^ shows these corrections as a 
functi on of luminosity, which agree broadly with the values 
in e.g. Richards et al. (2006c) over the luminosity range they 
consider. 

We have generally followed Ma rconi et alJ d2004t> in cal- 
culating this spectrum, with a more detailed treatment of the 
optical/IR and a more recent determination of a ox , but this 
level of detail is ultimately not required. The critical de- 
pendence is that of ctn X on luminosity. Adopting the median 
IRichards et alJ ll2006d) spectrum and rescaling the spectrum 
upwards of 0.1 keV according to the observed a ox for dif- 
ferent L„(2500 A) and corresponding bolometric luminosities 
yields a similar model spectrum and essentially identical con- 
clusions. Likewise, rescaling our median spectrum according 
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to Equation (21) in iMarconi et a Il ll200l gives a similar re- 
sult. Ignoring the dependence of bolometric corrections on 
luminosity altogether, however, is problematic. In the op- 
tical/UV, it is not as serious, as the bolometric correction 
changes by only ~ 20% from Mb = -23 to -30. However, 
over a ~ 4dex interval in L2-iokeV comparable to the usual 
observed baselines, the typical hard and soft X-ray bolomet- 
ric corrections change by more than an order of magnitude 
(primarily as given by the a ox -luminosity relation). 

Finally, even accounting for the dependence of spectral 
shape on luminosity, objects with a given luminosity do not 
all have identical spectra and bolometric corrections. It is im- 
portant to account for the dispersion in spectral shapes at a 
given L. In general, there will be two components, a corre- 
lated and uncorrelated dispersion. If, for example, there is a 
different spectral slope or value of a ox , then the bolometric 
correction in certain bands will be larger, but the correction 
in other bands must be smaller. Also, because the optical/UV 
contributes a larger fraction of the bolometric luminosity than 
the X-ray, the resulting dispersion in the B-band bolometric 
correction from different values of a ox will be smaller than 
the resulting dispersion in the X-ray bolometric corrections. 

We estimate these dispersions in the power-law compo- 
nents of our modeling from the observed distributions, as- 
sumin g that they ar e normally distributed; a ao « 0.125 
( Rich ard's et alJ 120031 compa ring the mean an in each ob- 
served quarti le), or = 0-30 (L Tozzi et aQ 120061) . and a am = 
0.075-0.14 dSteffen et alJl2006l) . Strictly speaking, a larger 
bolometric correction in one band implies a smaller integrated 
correction in others, but there are sources of scatter which 
introduce this effect weakly, such as e.g. variations in line 
strength in a given band or observational uncertainties in the 
luminosity in the band. Fitti ng to t he distribution of bolo- 
metric corrections in iRichards et al.l ll2006cl) . after account- 
ing for the luminosity distribution of the sources and corre- 
lated dispersions above yields a best-fit uncorrelated disper- 
sion component ~ 0.1 dex, consistent with observational esti- 
mates of th e scatter in band luminosities owing to these effects 
fe.g..lErvis et all 1 9941 iGeorge et alJl998HVanden Berk et al. l 
1200 1 | iPerola et al.ll2002l IRichards et al.ll2003t iHopkins et al.1 
2004). By fitting to a number of Monte Carlo realizations of 
the spectra as a function of luminosity given these dispersions, 
we can quantify the effective dispersion in the bolometric cor- 
rection in different bands as a function of luminosity 

0io g (£/£boi) = o"i (^boi/10 9 L Q )^+cr2 (3) 

with (ai,(3,a 2 ) = (0.08,-0.25,0.060) in the fi-band, 
(0.07, -0.17, 0.086) in the IR (15 /xm), (0.046, 0.10, 0.080) in 
the soft X-ray, and (0.06, 0.10, 0.08) in the hard X-ray. Here, 
a 1 and (3 roughly describe the correlated component of the 
dispersion, 02 the uncorrelated component. For typical bright 
quasars (Lboi ~ 10 13 L Q ), this reflects the fact that a significant 
(> 5%) fraction of quasars in any band can have their bolo- 
metric luminosities mis-estimated by a factor ~ 2 or more by 
a simple bolometric correction (even one that accounts for the 
luminosity-dependent spectral shape). These dispersions as a 
function of luminosity are plotted in Figure[T] 

2.2. The Observed Column Density Distribution 

In order to convert an observed luminosity function to a 
bolometric luminosity function, we must correct for extinc- 
tion in the different observed bands, which requires the adop- 
tion of an observed column density distribution. Essentially, 
the probability of observing a quasar of a given bolometric 



luminosity at some observed luminosity in a given band must 
account for the probability of extinction or attenuation. 

We consider three cases. First, our fiducial model adopts 
the luminosity-dependent observed column density distribu- 
tion fr om the hard and so ft X-ray observatio ns of lUeda et al.l 
(2003). We also follow lUeda et alJ ll2003l) and include an 
equal fraction of Compton-thick objects with Na > 10 24 cm" 2 
to that with Nu = 10 23 - 10 24 cm . The evidence for this in 
lUeda et ail (120031) is tentative, but it produces good agree- 
ment with the distribution of Comp ton- t hick c olumn densities 
subseq uentl y reported by Treis ter et al.l ( 120041) . iMainieri et alJ 
(2005), and iTozzi et alJ (|2006) and is consistent with upper 
li mits to the obscured fra ction from the mid-IR observations 
of Richards et alJ (l2006ci) . Recent very hard X-ray and soft 
gamma-ray (~ 20-200keV) Swift BAT and INTEGRAL ob- 
servations of local AGN, sensitive even to Compton-thick 
sources, confirm both a similar Compton-thick fraction and 
dependence on luminosity (demonstr ating also that this trend 
does not owe to selection effects) (iMarkw^ rdt et al] 120051 
iBeckmann et alJ l2006albl iBassani et alJ 120061 but see also 
Wang & Jiang 2006). Note that this yields a maximum 
Compton-thick fraction of ~ 30% (correcting the observed 
number density in a sample not sensitive to Compton-thick 
objects by a factor of 1 .4), at low luminosity. This does not 
imply a uniform factor of 1 .4 correction to the quasar space 
density, as the Compton-thick fraction depends on luminosity 
in the same manner as the entire column density distribution, 
and the fraction of Compton-thick sources at high luminosity 
will in general be much smaller. 

Alternatively, we consider a constant (luminosity- 
independe nt) column de nsity distribution, again adopting that 
from lUeda et alJ (12003) (fitted to their observations assum- 
ing no luminosity dependence). As the opposite extreme, 
we employ the colum n density distribution determined in 
lLa Franca et al.l ( 120051) . which depends on bot h luminosity 
and redshift. We discuss these models in § 12.31 but find that 
they are unable to produce a self-consistent set of luminosity 
functions in the different observed bands. 

Given an A^h distribution, we calculate the extinction at X- 
ray freq uencies using the photoel e ctric absorption cross sec- 
tions of Morrison & McCammon ( 1983) and non-relativistic 
Compton scattering cross sections. In the optical and mid- 
IR, we adopt a canonical gas-to-dust ratio (Ab/Nh)mw = 
8.47 x 10~ 22 cm 2 and Small Magellanic Cloud-like reddening 
curve from Pei ( 1992), which observations su ggest is appro- 
priate for the majority of r eddened quasars (Richards et alJ 
120031 IHopkins et all 120041 Ellison et aQ 120051) . although a 
Milky-Way like reddening curve changes the optical depth by 
only ~ 5 - 10% at the wavelengths of interest (excluding the 
2100 A "bump"). 

Given a bolometric QLF and the observed column density 
distribution, we can then convolve over the distribution of col- 
umn densities and spectral shapes at each bolometric lumi- 
nosity L to infer the implied distribution of luminosities that 
should be observed in a given band. In other words, knowing 
the probability of some intrinsic spectral shape and interven- 
ing column density, we determine the probability of observing 
quasars with an intrinsic L at some observed luminosity in the 
observed band(s). 

It is a surprisingly good approximation to this full convo- 
lution to adopt an "obscured fraction" as a function of lu- 
minosity; i.e. converting the luminosities of the bolometric 
QLF to luminosities in the observed band using the appro- 
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priate bolometric corrections and then multiplying by some 
"observable fraction" f(L) to correct for the effects of extinc- 
tion (and the dispersion in spectral shapes). Fitting to this 
f(L) from our full modeling yields a useful function in com- 
paring optical, soft X-ray, and hard X-ray observations; more 
directly applicable than the typically calculated fraction with 
Nu > 10 cm" 2 . This can be conveniently parameterized as a 
power law, 



AD 



<KL[L. 



lO^ergs" 1 



(4) 



where L,- is the luminosity in some band, and L is the cor- 
responding bolometric luminosity given by the bolometric 
corrections of Equation (0. This gives values (f^, (3) of 
(0.260,0.082) for L, = L B (4400 A), (0.438,0.068) for L, = 
L, R (15 /xm), (0.609,0.063) for U = L sx (0.5-2keV), and 
(1.243, 0.066) for L-, = L HX (2- lOkeV). Note that this "ob- 
servable fraction" can exceed unity, because the scatter and 
luminosity dependence of the bolometric corrections signifi- 
cantly changes the shape of the bright-end QLF in the X-ray 
bands (see Figure 0] below). These simple rescalings are ro- 
bust for the B and IR bands, with weak dependence on the 
shape of the bolometric QLF and dispersion in bolometric 
corrections, but the effects above make the observable frac- 
tion in the X-ray bands sensitive to the shape of the QLF; if 
greater accuracy is required, a more robust fit across a variety 
of QLF shapes gives j3 = 0.03572, 0.03472 for Lsx and Lhx, 
res pectiv ely, for an arbitrary bright-end QLF slope 72 (defined 
in § 13. 11 1. We caution that these are crude approximations, but 
the above equations can be used for rough conversion of ob- 
served QLFs to bolometric QLFs and vice versa, and for av- 
erage conversions between bands (multiplying any two such 
conversions together appropriately). 

2.3. Comparison of the Bolometric QLF Under Different 
Assumptions 

Given some model for the quasar spectrum and column 
density distribution, we can calculate the bolometric QLF 
from observations in some bands. In detail, we convolve a 
given bolometric QLF over the distribution of intrinsic spec- 
tral shapes and column densities, which yields the expected 
luminosity distribution in some observed band. By integrat- 
ing, if necessary, over the appropriate redshift and luminos- 
ity intervals, we can directly compare this to the binned QLF 
measurements at each frequency, luminosity, and redshift. We 
minimize the x 2 °f m is estimate in relation to the observed 
QLF in all bands, 
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to determine the best-fit bolometric QLF. Before applying this 
broadly, we would to test our description of the observed spec- 
tral shape and column density distributions. 

In Figure |2 we show (left panels) the bolometric QLF de- 
termined in this manner at several redshifts from our full mod- 
eling. At each redshift, we compare with QLF observations 
from Table[0with overlapping redshift intervals. We plot the 
binned QLF measurements rescaled to bolometric points with 
respect to the best-fit double power law bolometric QLF at 
that redshift (see §|3]below). In other words, convolving our 
best fit bolometric QLF with e.g. the quasar spectrum and col- 
umn density distribution predicts a number density of quasars 
«mdi for a given observed bin. Comparing this to the actual 



number observed, n b s , fixes the ratio « bs/«mdi- In Figure [2] 
we show the agreement with all bands simultaneously by plot- 
ting n bs/«mdi times the best-fit bolometric QLF at each lumi- 
nosity as the colored points. Given our full SED and obscura- 
tion model, the inferred bolometric QLF from each observed 
waveband agrees well. We quantify this directly, showing in 
each panel the x 2 l v statistic corresponding to the probabil- 
ity that all of the observations from the different bands derive 
from a single (technically double power-law, although this as- 
sumption only weakly changes the x 2 1 v ) bolometric QLF. 

If, however, we simplify our modeling of observed quasar 
spectra or column densities, the consistency between the QLF 
implied from the different wavebands is broken. First, we 
consider the prediction if all quasars ha d a single s pectra l 
shape by adopting the model spectrum of E lvis et all \ 19941) . 
with no dependence on luminosity. The optical, IR, and X- 
ray observations are no longer consistent with one another, 
and the x 2 l v at eac h redshift shown rises by a factor ~ 4 
to an unacceptably high value. Hard and soft X-ray mea- 
surements are still consistent, as they are in both cases re- 
lated by a relatively straightforward power law. These gen- 
eral conclusions are unchanged regardle ss of th e exact quasar 
spectrum adopted. In Richards et al. (2006c), mean spec- 
tra are computed for the entire quasar sample as well as for 
several sub-samples: the most optically luminous/dim, opti- 
cally red/blue, and IR luminous/dim halves (divided at the 
median values log[L opt /(ergs -1 )] = 46.02, A(g-i) = -0.04, 
and log[LiR/(ergs -1 )] =46.04). Figure |3]compares the bolo- 
metric QLF as in Figure |2] at z = 1. For clarity, we show 
only the best-fit double power- law rather than the b inned 
points. Consideri ng the different Richards et al. (20063) an d 
lElvis et all ( 119941) mean spectra, no single quasar spectrum 
yields a consistent bolometric QLF at all luminosities (i.e. 
there is no "effective mean" spectral shape). Taking the mean 
spectrum of the most optically or IR bright quasars, unsur- 
prisingly, yields a consistent bolometric QLF at the bright 
end (above the break). The optically and IR dim spectra, on 
the other hand, are appropriate for lower-luminosity quasars 
(those near and just below the break), and produce consis- 
tent bolometric QLFs in this regime. Note that this does 
not extend to th e lowest luminosities s hown, as the "optically 
dim" objects in Richards et al. (2006c) are still much brighter 
(Lboi > 10 45 ergs _1 rs 3 x 10 11 Lq) than the lowest luminosi- 
ties plotted in Figure[3]and probed by X-ray samples. 

We next examine different assumptions for the column den- 
sity distribution. In Figure|5] we adopt the full spectral model 
but consider the lum inosity-in dependent mean column den- 
sity distribution from lUeda et al.l (12003): i.e. a constant ratio 
of obscured to unobscured AGN ~ 2 : 1 , similar t o that found 
locally (e.g- lRisaliti et al.lll999HHao et al.1l2005l) . Again, the 
agreement is broken: qu asars at low l umino sities, where the 
number of quasars in the lUeda et alJ (12003) sample is large, 
dominate the mean A^h distribution, and as a result this con- 
stant obscured fraction provides an acceptable (although still 
less good) fit to the low-luminosity data. However, at high 
luminosities, the results derived from different wavebands di- 
verge. The total x 2 l v unambiguously rises (factor ~ 2) at all 
redshifts, with an increasing discrepancy at higher luminosi- 
ties (factor ~ 3 increase in x 2 l v above the QLF break). The 
disagreement between different bands changes significantly 
as a function of luminosity, so re-normalizing the constant 
obscured fraction cannot make the luminosity functions con- 
sistent. 
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FIG. 2. — The bolometric QLF from the observations in Table HI at various wavebands (optical: green, soft X-ray: blue, hard X-ray: red, IR: cyan; see 
Table HI for the plotting symbols). For ease of comparison, the points in the different bands are rescaled as described in the text to indicate the bolometric 
QLF implied by each set of measurements. Left panels show the bolometric QLF adopting our full (luminosity-dependent) bolometric corrections and column 
density distributions (see also Ueda et al. 2003; Marconi et al. 2004; La Franca et al. 2005), at several redshifts (as labeled). Center-left shows the result with 
the full column distribution, but adopting the constant bolometric corrections of lElvis et alj 119941) . Center-right uses the full bolometric corrections, but a 
luminosity-independent column density distribution (obscured fraction). Right uses the full bolometric corrections and a strongly redshift-dependent column 
density distribution from La Franca et al. 1 2005). Each panel shows the reduced \ 2 f° r the assumption that the data yield a consistent bolometric QLF at each 
redshift. The latter three assumptions do not yield a consistent bolometric QLF at each redshift, unlike for the observationally derived luminosity-dependent 
bolometric corrections and column density distributions. 



Finally, we adopt a column density distribution which 
evolves strongly with redshift (and luminosity ). Specifically, 
we adopt the fit from lLa Franca et alJ (120051) to the Nn dis- 
tribution with maximal redshift evolution, similar to the red- 
shift evolut i on im plied by the column density distributions in 
iTozzi et ID (|2Q06) (fitted assuming maximal redshift depen- 
dence and no luminosity dependence), and similar to the red- 
shift evolution of obscured fractions implied by some X-ray 
background synthesis models (e.g., Comastri et al. 1995; Gilli 
et al. 1999; but see also Treister & Urry 2005). The resulting 
bolometric QLF is shown in the right panels of Figure|2] 

This fit is reasonable at low and moderate redshifts z < 1 
(unsurprisingly, where the observed samples are best con- 



strained). Extrapolated to higher redshifts, the consistency 
between different bands breaks down, increasing \ 2 l v by an 
unacceptable factor ~ 3 by z > 2. Therefore, while our anal- 
ysis implies that the obscured fraction must not be constant, 
it is also clear that any variation is primarily a luminosity, 
and not a redshift dependence. However, more moderate red- 
sh ift evolution such as tha t suggested bv lBarger et alJ J2005I) 
or lBallantvne et all 1^2006), with obscured fractions evolving 
by only ~ 20% from z = — 1 (with no further evolution at 
higher redshifts), is consistent with our analysis, and will have 
only a small impact on the cumulative QLF evolution. Such 
evolution, being a small effect, would be better constrained 
by direct measurements of the column density distribution at 
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FIG. 3. — Comparison of the best-fit bolometric QLFs to the observations in Figure[2]at z ~ 1, fitted independently to the optical (green), soft X-ray (blue), 
hard X-ray (red), and IR (cyan) data sets. Upper left panel shows the result using our full (luminosity-dependent) bolometric corrections, with the black line 
the bolometric QLF from fitting to the data at all wavelengths simultaneously (reproduced as the dashed black line in the other panels), subsequent panels adopt 
constant (luminosity-independent) corrections from Richards et al. 1 2006c) based on the mean sp ectrum of th e complete (all), blue, red, optically bright, optically 
dim, IR bright, and IR dim sub-samples considered therein. Lower right panel adopts the Elvis et al. 1 1994) mean spectrum. Accounting for the dependence of 
spectral shape on luminosity yields a consistent bolometric QLF across all frequencies and luminosities, whereas considering a single mean spectrum is generally 
appropriate only for narrow luminosity and color intervals. 



moderate redshifts. 

Likewise, although not shown in Figure [2] we find no evi- 
dence for any dependence of bolometric corrections on red- 
shift, in agreement wi th a number of d i rect measurements 
of quasar SEDs ("e.g.. Elv is et alJ [ 19941: fVa nden Be rk et all 
1 20011 iTelfer et alJl2002l IVignali et alj|2003l: iFan et alJl2003t 
Steff en et alJ 120061 iRichards et alJ 12006c ; Shemmer et al. 
2006). If the variation of e.g. a ox were primarily a redshift 
as opposed to luminosity dependence, or if there were e.g. a 
strong depende nce of X-ray photo n index on red shift as sug- 
gested bv iKuhn et ail J2001I) and lBechtold et alJ <2003l) . then 
our inferred bolometric QLFs from different bands would di- 
verge with redshift. Instead, they appear to be self-consistent 
given the redshift-independent spectral model at all z = 0-6. 
Given the limitations of present data, however, it is still possi- 
ble that less well-constrained portions of the spectrum evolve 
differently. For example, the IR SED could, in principle, 
evolve differentl y with luminosity or redshift from the opti- 
cal SED (but see lJiang et al|2006bl who find simi lar near-IR 
SEDs at high redshift), and Maiolino et al. (2004) have sug- 
gested that extinction curves may evolve at z > 4-5 (although 
their proposed extinction curves change the optical depth at 
the wavelengths of interest by < 10%). There is no strong 
evidence for such additional evolution in the samples we con- 
sider, but these caveats should be considered in any extrapo- 
lation of our fitting to less well-sampled frequencies and red- 
shifts. 



2.4. The Relation of Bolometric to Observed QLFs 

We briefly examine the relation between the shape of the 
bolometric QLF and the observed QLF at different frequen- 
cies. Figure^ shows the fit to the bolometric QLF from Fig- 
ure 13 at z ~ 1, and the resulting observed QLF in different 
bands. Although our fitting considers all observations in the 
bands given in Tabled we have renormalized all optical ob- 
servations to the B-band for plotting purposes (likewise for 
the other bands). We plot all observations in the same units 
(B-band and 15 [im show vL v ) for the sake of direct compari- 
son. 

The optical/UV spectrum dominates the bolometric lumi- 
nosity (at least in bright, unobscured quasars), and therefore 
the optical bolometric correction depends only weakly on lu- 
minosity. Consequently, the bright-end slope of the optical 
QLF is essentially identical to that of the bolometric QLF, off- 
set by a nearly constant (when calibrated for optically bright 
quasars) bolometric correction. However, at the faint end, ob- 
scuration significantly flattens the optical QLF. In the hard 
X-ray, this is reversed. The faint-end slope is relatively un- 
affected by obscuration and directly probes the bolometric 
faint-end, but the change in bolometric corrections (as well 
as the scatter in bolometric corrections at fixed L) alters the 
bright-end shape considerably. The IR, unsurprisingly, pro- 
vides in many ways the most direct estimate of the bolometric 
luminosity distribution, but is the most poorly constrained by 
current observations. The soft X-ray, while yielding some of 
the tightest observational constraints and critical for evaluat- 
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ing e.g. the UV contributions of quasars, is affected at the 
faint and bright ends by obscuration and changing bolometric 
corrections, respectively. 

2.5. The Contribution of Quasar Host Galaxy Light 

We have so far ignored the contribution of quasar host 
galaxies to their observed luminosities. Observed luminos- 
ity functions, however, do not generally remove this host 
contamination. In the X-ray bands, this is not of course 
expected to be a serious ^concern, and in the emission-line 
luminosity functions of Hao et al. (2005) we consider their 
more conservative AGN cut which should eliminate much of 
the contamination from star-forming nuclear regions in the 
hosts. However, this contribution may still be problematic 
at optical and IR wavelengths. If we knew the black hole 
masses of quasars in the observed QLFs, it would be straight- 
forward to at least estimate an average host contribution or 
bias, as there is a well-established correlation between black 
hole mass and host galaxy op tical luminosity (or mass) (e.g., 
iKormendv & Richstone 1995). Lacking this information, we 
must adopt some estimate of the observed quasar Eddington 
ratio (L/Lndd) distribution to convert from a bolometric lumi- 
nosity to black hole mass and corresponding host luminosity. 

In Figure |5] we consider several simple, representative 
cases to estimate the possible effects of our neglecting host lu- 
minosities. We compare the Z?-band QLF neglecting host light 
(our fiducial model) and alternatively, assuming all quasars 
are radiating at the Eddington luminosity. This gives a black 
hole mass for each L, and we use the b lack hole mass- 
host fl -band luminosity relationship from Mar coni & Hunt! 
fc00 3 |) to determine the corre sponding host contribution (see 
alsoEanden Be rk et alJ l2006. who measure a similar L/Lndd- 
dependent relation directly in AGN). The quasar luminosi- 
ties are corrected and attenuated as before, and then the (un- 
attenuated) host luminosity is added. 

Strictly speaking, this is not necessarily self-consistent, as 
the observational calibrations of bolometric corrections do not 
necessarily remove host galaxy contributions. However, our 
construction of the intrinsic quasar SED should effecti vely ex- 
clude host contributions (compare e.g. Richa rds' et alj f2006c) 
who explicitly remove a rough lower limit to the host contri- 
bution), so long as there is not a substantial host contribution 
at ~ 2500 A which would significantly affect the determina- 
tion of a ox (altho ugh this is not expect ed at these frequencies, 
and the relation in St effen et al.N200 6) is calibrated from only 
Type 1 AGN). In any case, these effects are second-order to 
those shown in Figure|5] 

We also consider the case where all quasars radiate at an 
Eddington ratio L/Lem = 0.3 and L/Lgdd = 0.1. As an alterna- 
tive simplification, we consider the case if the QLF is a pure 
Eddington ratio sequence with L/Lem = LjL^_ (L* is the fitted 
QLF break luminosity). It makes no difference if we allow 
L/Lem > 1 or fix the maximum L/Lndd = 1, as the host contri- 
bution is minimal in either case. We show this for the fi-band, 
but the effects in the mid -IR are similar (given e.g. the typi - 
cal early-type spectrum of lFioc & Rocca-VoImerangel l ll997l) : 
effects are also similar for other host types). 

Figure|5]demonstrates that including quasar host galaxy lu- 
minosities should be a small effect, except at the lowest lu- 
minosities Mb > -23 (bolometric L < 10 I2 L Q ). This is un- 
surprising, as at bright observed Mb, obscuration is small 
and Eddington ratios are expected to be high. In general, 
this implies a negligible impact on our subsequent calcula- 
tions, as the deepest optical observations only just probe the 



range in which host luminosity contamination becomes sig- 
nificant (the larger baseline at low luminosities comes from 
hard X-ray observations, which do observe significant num- 
bers of optically normal galaxies hosting low-luminosity X- 
ray AGN, e.g. Barger et al. 2005). At higher redshift, flux lim- 
ited samples are further removed from these luminosities, and 
Eddington ratios are expected to be unifor mly high (at leas t 
in optically-selected samples, e.g. Mc Lure & Dunlopll2004l) . 
Moreover, observations of the relation between black hole 
mass and host luminosity at high redshift, albeit considerably 
uncertain, suggest that host luminositi es either remain con- 
stant or decrease at a given Mbh (e.g. JPeng et al.l 2006). 

However, this strongly cautions the interpretation of deeper 
optical, near- or mid-IR surveys, especially at low redshifts 
z < 0.5, where Eddington ratios may well be low (< 0.1). 
If such surveys seek to probe luminosities < 10 12 L Q , care- 
ful consideration of the contribution from quasar hosts is nec- 
essary to determine the actual quasar contribution at low lu- 
minosities. For example, Figure |5] demonstrates that the ob- 
served faint-end QLF slope in these bands can be significantly 
biased by host contamination at magnitudes lower than those 
currently probed, in a manner sensitive to the Eddington ra- 
tio or black hole mass distribution. These effects will be even 
more pronounced in the near-IR, as the ratio of host to quasar 
luminosity in unobscured object s has a typical maxim um at 
v ~ 1 .6/im (see e.g. Figure 1 1 of Richards et al. 2006c). 

3. THE BOLOMETRIC QLF 
3.1. The Bolometric QLF at Specific Redshifts 

We combine the binned QLF measurements to examine the 
QLF "at" a single redshift. Of course, the observations are 
not made over identical redshift intervals, so we re-normalize 
each to the same redshift (as « bs/imdi) with the best-fit ana- 
lytic QLF fit from the same (or appropriate companion) paper. 
This introduces some additional model dependence, so these 
fits should be considered as heuristic, but they will inform our 
subsequent choice of functional forms for more properly fit- 
ting to the redshift dependence of the bolometric QLF. 

We follow standard practice and fit the QLF to a double 
power law 
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with normalization break luminosity L», faint-end slope 
71, and bright-end slope 72. Note that conventions for the 
double power law are often different in optical and X-ray anal- 
yses; for optical QLF s typically the double power law is de- 
fined in terms of (e.g., Peterso rll997tlCroom et all2004l) 
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dL (L/L*)- a + (L/L*)-P' 
or per unit absolute magnitude 

d$ = & 

AM l00Ma+l)(M-M t ) + IQ0A(J3+1)(M-M,) ' 

which in our notation gives a = -(71 + 1), j3 = -(72+ 1), </>'* = 
<£,/lnl0, and $f = O.40*. 

The resulting bolometric QLF at several redshifts is shown 
in Figures [6] & At each redshift, all binned observations 
from Tableau overlapping redshift intervals are shown. The 
QLFs from optical, soft and hard X-ray, IR, and emission-line 
measurements yield a self-consistent bolometric QLF at all lu- 
minosities and redshifts probed, when a full modeling of the 
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FIG. 4. — The best-fit bolometric QLF at z ~ 1 (black lines), with the resulting observed QLF in each of several bands: optical (green), soft X-ray (blue), hard 
X-ray (red), and IR (cyan). Each corresponding panel plots the compiled observations in TablelTlfor the appropriate redshift and frequency (points). Symbols 
denote the parent sample as given in Table fTI The observations at each band are consistently produced from a single bolometric QLF, and together provide 
strong constraints on the QLF shape. Different wavebands accurately represent different aspects of the bolometric QLF shape. With only a weakly L-dependent 
bolometric correction but significant effects of obscuration, the optical QLF faithfully traces the bolometric bright end shape, but is flatter at faint L. Hard X-rays, 
conversely, are weakly affected by obscuration but can have a strongly luminosity-dependent bolometric correction, reproducing the faint end shape but being 
steeper at bright L. The IR QLF better follows the bolometric shape, but is currently least well-constrained. 



observed spectral and Nn distributions is employed. Further- 
more, the combination of this number of observations pro- 
vides a large luminosity and redshift baseline, sampling the 
the faint end, break, and bright-end slope to redshifts z ~ 4-5. 

We plot the best-fit parameters from these fits at a number 
of redshifts in Figure[8](errorbars show formal lcr uncertain- 
ties from the fits), and list them at several redshifts in Table[2] 
The normalization 0* is roughly constant, while the break lu- 
minosity (which is quite tightly constrained for most of the 
redshift range) evolves by ~ 2 orders of magnitude. Note that 
the points where 0* appears to deviate from being constant 
also have discrepant values — the degeneracy between the 
two is such that the value of 0* is consistent with being con- 
stant at all z for smooth evolution in L* . There is an indication 
of evolution in both the faint-end and bright-end slopes, which 
we discuss below. The integrated bolometric luminosity den- 
sity is well constrained, with the largest uncertainties only 
~ 0.15dex. We determine the luminosity density by integrat- 
ing the best-fit luminosity function to L = 0. At most redshifts 
z > 0.5 the faint-end slope is relatively shallow, and choosing 
instead a cutoff at e.g. L = 10 8 - 1O 9 L0 changes the integrated 
luminosity density by < 10%. At the lowest redshifts, how- 
ever, the faint-end slope is steep, and at z = where this is 



most pronounced, the luminosity density is ~ 15% lower if 
we truncate at L = 10 8 L Q (~ 25% lower for L = 1O 9 L ). This 
sensitivity to the steep faint-end slope at low-z is the reason 
for the relatively large uncertainty in the luminosity density 
at low redshift. At high-z the uncertainty owes to the limited 
amount of data. In any case there is a well-defined peak in the 
luminosity density at z = 2. 154 ± 0.052 (formal error from fit; 
we expect a systematic error ±0.15 from choices in sampling 
and binning the observations), well outside the range where 
either of these systematic concerns is problematic. 

3.2. Analytic Fits as a Function of Redshift 

We characterize the QLF as a function of redshift by adopt- 
ing a standard pure luminosity evolution (PLE) model, where 
the bolometric QLF is a double power law at all z, with con- 
stant 71 , 72, and , but an evolving L* . We allow L. t to evolve 
as a cubic polynomial in redshift, 

log U = (log )o + k L . , £ + k L , 2 e 2 + k L . 3 £ 3 , (9) 

where 




Here, ki\, ki2, and k^3 are free parameters, and we set 
z re f = 2 (which roughly minimizes the covariance in the fit). 
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FIG. 5. — The observed B-band QLF (points) at z ~ 1, and that determined 
from the best-fit bolometric QLF in Figure|4](solid line) ignoring contribu- 
tions to the observed luminosity from quasar host galaxies. Also shown is the 
fi-band QLF determined from the same bolometric QLF, but including the ex- 
pected host galaxy contribution to the observed luminosities (normalized for 
a given Eddington ratio with the observed relation between black hole mass 
and host galaxy B-band luminosity from Marconi & Hunt 1 2003)), assuming 
all quasars have constant L/LeiM = 1-0, 0.3, 0.1 (dot-dashed, long-dashed, 
dotted lines, respectively), or assuming L/Leu =L/L„ (short dashed line; L» 
is the QLF break luminosity). Host galaxy contributions should be negligible 
at most luminosities and redshifts for the observed bands we consider, but 
will be important for future optical, near- and mid-IR surveys which probe 
fainter luminosities L < 1O 12 L (M B > -23). 



The cubic term is demanded by the data (A% 2 <~ 600 on its 
addition), but higher order terms in £ are not (A\ 2 < 1). Since 
this model includes the evolution with redshift, we can simul- 
taneously fit to all the data sets in Tabled each over the ap- 
propriate redshift intervals of the observed samples. 

The best-fit PLE model parameters are given in Tableland 
plotted as a function of redshift in Figure|S] and the resulting 
QLF is shown in Figures[6]&Gl Although this provides a rea- 
sonable lowest-order approximation to the data, it fails at the 
faint end, underpredicting the abundance of low-luminosity 
sources at z < 0.3 and overpredicting it at z > 2, and the fit is 
poor at z > 5 with much too steep a bright end slope. Over 
the entire data set, the fit is unacceptable, with \ 2 ~ 1924 for 
v = 510 degrees of freedom. A pure density evolution (PDE) 
model fares even worse, with X 2 /v = 3255/510 (although for 
completeness we provide the best-fit PDE parameters in Ta- 
ble 0}, unsurprising given that nearly every observed data set 
which resolves the break in the QLF favors the PLE form 
(e.g.. iBovle et alJl2000l iMivaii et alj|2000l lUeda et alJl2003t 
ICroometalJ2 004: Richards J etall2005l) . 

As discussed in § ^ many recent studies have found ev- 
idence for evolution in the faint end beyond that predicted 
by PLE, with the density of lower-luminosity sources peak- 
ing at lower redshift than the density of higher-luminosity 
sources, and have fit this trend with a luminosity-dependent 
density evolution (LDDE) model (Schmidt & G reed fl983l) . 
while high-redshift samples have suggested evolution in the 
bright-end QLF shape, confirmed robustly as a flattening in 
7? for the first time in a homogeneous sample by the SDSS 
jRichards et alJ2006d) . We wish to incorporate this additional 
evolution into our model. 

For comparison with these results and future X-ray surveys, 
we fit to an LDDE form allowing maximal flexibility of the 
parameters, 



{L/L*Vi+{L/U)n 
The density function e& is given by 



e&(L, z). 



(Z < Z C ) 



with 



and 



(l+z) pl 

(l+z c )" l [(l+z)/(l + z c )Y 2 (z>Zc) 
z c (L) = 



z c ,o(L/Lc) a (L<L C ) 



Z c ,0 



(L > L c ) 



(12) 



(13) 



(14) 



/? l(L) = /9l 4 6 + A[log(L/10 46 ergs- 1 )] (15) 
P 2(L) = p2 46 + fa [log(L/ 1 46 erg s" 1 )] . (16) 

Note that some authors ("e.g.. lHasinger et afll2005h adopt an 
alternative normalization convention in terms of A44 = 4>(L = 
10 44 ergs _1 , z = z c ) and z c ,44 = Zc(L = 10 44 ergs _1 ), but these 
choices are not particularly convenient for the bolometric 
QLF. The 11 free parameters in this fit are then L*, <p^, 71, 
72, z c ,Q, L c , a, pl46, p2^, /?!, /?2- Their best-fit values are 
given in Table|4] 

The LDDE form can effectively describe evolution in the 
faint-end QLF slope, but it does not allow for evolution in the 
bright-end slope. It also has a tendency to introduce a "second 
break" in the faint end of the QLF (i.e. if the faint end flattens 
there is often some L below which it rises steeply again), for 
which we see no evidence. Ultimately, the improvement over 
the PLE fit is highly significant, with x 2 /v = 1389/507, but 
there is still room for substantial improvement. 

Therefore, we instead consider the PLE form above, but 
allow both the bright- and faint-end slopes to evolve with red- 
shift. For the faint-end slope, using the fitted points at each 
z to inform our choice of functional form, we model 71 as a 
power-law in redshift, 



7i =(7i)oxlO^> c 
= (7i)o(t^)^' 



(17) 
(18) 



where again z re f = 2 is fixed. Allowing this dependence (while 
still holding 72 constant) significantly improves the quality of 
the fit relative to the PLE model, reducing \ 2 by &\ 2 ~ 500 
with the addition of one parameter, to \ 2 l v = 1422/510. The 
values for this fit are given in Table [3] There is no evidence 
for higher-order terms ( A\ 2 < 1 for the addition of a second- 
order power in £). We have also tested different functional 
forms, and do not find any which provide a significantly better 
fit. This choice also has the advantage that it extrapolates to a 
flat 71 — > (as opposed to a negative, likely unphysical) slope 
at high redshift. 

Parameterizing evolution in the bright-end slope is more 
difficult. There appears to be evidence for a steepening of 
the bright-end slope (increase in 72) from z ~ to z ~ 1.5, 
then a flattening with redshift. We find the best results for a 
double-power law of the form 



72 = (72)0 x 2 



2(72)0 



\ 1+Zref J \ l+Zrcf 



\k-i2-- 



(19) 



(20) 



<j>{L,z) = <j>{L, Q)e d (L,z) 



Here, £ 72 1 describes the rise of 72 with z at low redshift, and 
(11) 2 describes the fall of 72 with z at high redshift. The fac- 
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FIG. 6. — The best-fit bolometric QLF at each of several redshifts (left panels; shown as M m di/n bs), and the corresponding observed QLF in B-band (center-left; 
green), soft X-rays (0.5 -2 keV) (center; blue), hard X-rays (2 — lOkeV) (center- right; red), and mid-IR (15 /im) (right; cyan). Rather than add a series of panels 
for a single data set, the emission-line luminosity functions of Hao et al. 1 2005) are shown (orange) in the z — 0.1 hard X-ray panel (rescaled by « bs/«mdl: but 
equivalently directly converted to hard X-ray luminosities following Heckman et al. 2005). Lines show the best-fit evolving double power-law model to all 
redshifts (solid), the best-fit model at the given redshift (dashed), and the best-fit PLE model (dotted). Points shown are the compiled observations from TablefTI 
with the plotting symbols for each observed sample listed therein. 



tor of 2 is inserted such that if hy lt \ = £ 7 ,.2 = 0, the result- 
ing 72 = (72)0- The best-fit parameters, fixing 71 or allow- 
ing both slopes to evolve simultaneously, are given in Ta- 
ble [3] We find that k 12> \ is non-zero at high formal signifi- 
cance, but cannot say whether this steepening up to z ~ 1 .5 
is "real" in any robust sense; however, including this pa- 
rameter does significantly improve the accuracy of our fit- 
ting function in representing the binned data. The flatten- 
ing at high redshifts is more significant (~ la). Including 
these two parameters greatly improves the quality of the fit by 
A\ 2 ~ 600, giving x 2 /v = 1312/509 for a fit with constant 71 
or x 2 /^ = 1007 / 508 for a simultaneous fit including evolution 
in 71 and 72. Again, there is no evidence for terms describing 
further evolution. Although our final best fit does not cross 
this limit until the highest redshifts, we also formally enforce 
a lower limit 72 > 1 .3 (to prevent an unphysical divergence in 
the luminosity density). 

The best-fit parameters as a function of redshift (allowing 
all free parameters to vary) are plotted in Figure [S] There 
is a slight offset between the L* and 0, of this fit and that 
fitted at specific redshifts, but this owes to their covariance. 



This can be seen in the bolometric luminosity density, which 
accurately traces the fit predictions (note that the PLE model 
significantly overpredicts the luminosity density at z > 2, as 
it overpredicts the number of faint sources). Our restriction 
72 > 1 .3 is important in extrapolating beyond z ~ 5 to prevent 
a divergence in the number and luminosity density of bright 
sources. We also plot the predicted number density of quasars 
with Mb < -27 as a function of redshift, which agrees with the 
direct observations (expected, but nevertheless a reassuring 
consistency check). 

Figure|9]plots the number density of quasars integrated over 
various luminosity intervals, in various bands as a function of 
redshift, from the best-fit model and the compiled observa- 
tions in Table Although this information is contained in 
Figures |6] and Figure [9] nicely illustrates an essential trend 
captured by the LDDE or evolving double power-law forms, 
namely that the density of lower-luminosity sources peaks at 
lower redshift than that of higher-luminosity sources. The 
trend is evident in all bands we consider. 

We briefly consider an alternate model for the shape of the 
QLF, which illustrates several important points. Instead of 
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the double power-law of Equation (jSJl, we consider a modified 
Schechter function 

<IW) = 0* (L/L*)~ 71 exp { - (A) 72 }. (21) 

We have also, for example, adopted polynomials of arbitrar- 
ily high order (although the fits typically do not improve be- 
yond fourth order). In either case, the fit is quite similar at 
most luminosities (with similar \ 2 /v), implying that there is 
no dramatic shape dependence which we have not captured. 
However, such a fit does exhibit smoother curvature rather 
than a sharp break at L*, evi dence for which has been seen in 
some optical samples (e.g., Wolf et al. 2003; Richar ds et al J 
120051) . At the highest luminosities (> 1 dex above L*, typi- 
cally Lboi > \Q ia Lq) the implied number of quasars is an or- 
der of magnitude lower for these parameterizations than for 
the double power law prediction, and falls much more rapidly. 
The resulting observed luminosity functions are more sensi- 
tive to the estimated dispersion in bolometric corrections, but 
in either case the highest luminosity soft and hard X-ray ob- 
jects in Figures |6] & are substantially affected by quasars 
shifting into slightly larger bins of soft or hard X-ray lumi- 
nosity owing to different spectral shapes (e.g. the scatter in 
a ox ). It is important to account for this effect when attempt- 
ing to infer the number density of the most massive black 



holes and most luminous quasars, as a naive extrapolation of 
the median bolometric corrections applied to the most X-ray 
bright quasars implies extreme (and potentially unphysical) 
bolometric luminosities > 10 15 L Q (i.e. a > 3 x 10 I0 M© black 
hole at the Eddington rate). Multiwavelength observations 
of these particular objects and further study from large area 
surveys which do not have to bin in widely spaced luminos- 
ity intervals will be critical in breaking the degeneracies be- 
tween these fits to the intrinsic bolometric QLF and the double 
power-law form. 

We are unable to find any further dependences which sig- 
nificantly improve the best-fit QLF. Allowing the parameters 
describing the observed column density distribution and spec- 
tral shape, and their respective luminosity dependence, to vary 
simultaneously, yields only marginal improvement. It appears 
that the remaining scatter in the data does not mostly owe to 
a failure to capture some remaining dependence. In Tabled 
we list the x 2 l v for each sample with respect to the best-fit 
full model of the evolution of the QLF in Table|3] The agree- 
ment with most samples is good, and the largest, most well- 
constrained samples (with small typical <C 0.1 dex errorbars) 
give x 2 l v ^ 3, in each case comparable to or smaller than 
the reduced x 2 found by the respective authors in fits to those 
individual data sets. To the extent that the functional forms 
adopted by these authors cannot reduce the variance below 
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FIG . 8 . — The best-fit QLF double-power law parameters as a function of redshift. Points show the best-fit values to data at each redshift, dotted lines the best-fit 
PLE model, an d solid lines the best-fit full model (with cyan shaded range showing the 1<t uncertainty). Open diamonds in 72(2) show the bright-end slope fits 
from Richards et al. 1 2006b). Although PLE is appropriate for a lowest-order fit, both the bright and faint-end slopes evolve with redshift to high significance 
(> 6a ). Lower right shows the predicted number density of bright optical (Mb < —27) qu asars from the full fit (solid), compared to that observed in lCroom et all 
i 2004) (square), Richards et al. 1 2006b) (diamonds and dashed line), and Fan et al. 1 2004) (circle). 



this level, we would not expect to do better in the combined 
sample. 

There will also be some unavoidable variance introduced 
owing to systematic variations between independent data sets. 
Various observational calibrations, the model dependences in- 
herent in calculating a binned QLF, and most of all cosmic 
variance will all contribute to sample-to-sample differences. 
In fact, we find that allowing for even a small ^0.05 dex 
(10-15%) systematic normalization variance between sam- 
ples, most of the remaining scatter is accounted for, with a 
best-fit model improvement A% 2 ~ 500. We provide the val- 
ues of this fit in Table[5] but caution that we have increased the 
sample-to-sample variance by this amount uniformly, when in 
fact systematic effects such as cosmic variance will be smaller 
for large surveys such as the SDSS and 2dF than for the small 
deep fields of Chandra and XMM. Consequently, this under- 
weights some of the most well-constrained observations, and 
the fit results should be considered to be heuristic. 

4. THE BUILDUP OF THE BLACK HOLE POPULATION 

4.1. Model-Dependent Quantities 

We can gain further insight into the evolution of the 
QLF, albeit at the cost of some model dependence, by de- 
convolving the observed quasar luminosity function with 
a theoretical model for the quasar lig ht curve or lif etime. 
This method is well-e stablished (e.g . . ISalucci et alJ 119991 
| Yu & Tremaml 120021 iMarconi et al! 120041 iShankar et alJ 
l2004UYu & Lul2004l) . but most studies adopt simplified "toy" 
light curves with accretion at a constant (fitted) absolute rate 
or Eddington ratio and fitted lifetimes and duty cycles. We 



instead follow Hop kins et alJ d2006bl) . who derive physically 
motivated quasar light curves from simulations of merging 
galaxies that include black hole growth (e.g. Springel et al. 
2005b) and which are consistent with a large range of obser- 
vational constraints that cannot be reproduced by idealized 
models. This also removes the various fitting degeneracies 
- for a given bolometric QLF, the quasar light curves deter- 
mined in simulations yield a unique black hole mass function, 
cosmic X-ray background spectrum, and self-consistent black 
hole and host galaxy properties. 

Given a consistent model of the quasar lifetime/lightcurve, 
the observed bolometric QLF is given by the convolution of 
the rate of quasar formation or "triggering" with the differen- 
tial quasar lifetime, 



<KL) = / <KMbh) 



dt 



dlogL 



(L | M BH ) d \ogM BH (22) 



where 



4>{M BH ) = 



d$(M BH , t) 
d\ogM BH &t 



(23) 



is the rate of formation of black holes of a relic mass M B h 
at cosmic time t. Since dt /dlogL is completely determined 
in the simulations and analytical models of Hopki ns et alJ 
j2006albl) . we can fit to 4>(M BH ) in the same manner that we 
have fitted 4>{L). 

We assume a double power-law form for c/)(Mbh) at all z, 



(M Bn /M*)v<+(M Bll /M*)vi 



(24) 
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FIG. 9. — Total number density of quasars in various luminosity intervals (in log (Lband /[ergs - ]) as labeled) as a function of redshift, from the best-fit evolving 
double-power law model (lines) and the compiled observations in TablefT1(points of corresponding colors, symbols for each sample listed therein) in bolometric 
luminosity, 73-band, soft X-rays (0.5—2 keV), and hard X-rays (2— lOkeV). The trend that the density of lower-luminosity AGN peaks at lower redshift is manifest 
in all bands. 



with normalization <f>*, break M*, and faint-end and bright- 
end slopes rji and 772, respectively. We determine an analyti- 
cal fit to <P(Mbh, z) in the same manner as <f>{L, z). We allow 
logM* to vary as a cubic in £ just as L», and allow the high- 
mass slope 772 to vary just as 72, as well. There is no evidence 
for evolution of the low-mass slope 771 . In this formulation 
there is evidence for evolution in </>» above z ~ 2, so we fix it 
to a constant below z re f = 2, and allow it to evolve as a power- 
law above z re f, 



((f>*)o (Z < Zref) 

(0*)O[(1+Z)/(1+Zref)]** (z>Zref). 



(25) 



In Figure[H)] we plot the best-fit normalization <p* and break 
(characteristic mass in formation) M* as a function of red- 
shift. At high redshifts, objects build up rapidly, until z ~ 2, 
after which the rate of merger/black hole growth "events" flat- 
tens, and activity ceases in the most massive systems and 



rapidly moves to less massive objects at lower redshifts, per- 
haps driven by feedback mechanisms quenching activity in 
the higher-mass systems. Given the adopted lifetime mod- 
els, the observed faint-end <p(L) slope 71 is dominated by 
sources with L <C LEdd(^BH) and is determined by the quasar 
lifetime as a function of Mbh- Figure ^| shows the pre- 
dicted 71 (z) from the quasar lifetime model given the best-fit 
4>(Mbh) at each z, compared to the direct fits from Figure |8j 
the agreement is good despite rj\ being nearly constant and 
much flatter 771 fti 0.0- 0.2 at all redshifts (see also Figure 3 of 
iHopkins et aLH2006al) V 

4.2. Integrated Quantities 

If all black holes accrete with some constant radiative effi- 
ciency e r , 



Lbol = e,Mc z 



(26) 
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FIG. 10. — Best-fit parameters to </>(Mbh). the "birthrate" or formation rate of black holes at a given redshift, given the quasar lifetimes determined in the 
simulations of Hopkins et al. 1 2006a b). Top left shows d<f>/d\, the quasar activation / birthrate normalization. Top center shows M«, the break in the characteristic 
final black hole mass forming at a given z. Points are fitted at each redshift, solid lines show best-fit from </>(Mbh) as a function of redshift. Top right shows the 
faint-end slope of the observed QLF from Figure|8]at each redshift (black circles), with the predicted faint-end slope from the same quasar lifetime models (red 
crosses determined from the observations at each z). Bottom left shows the total black hole mass density of the Universe, from the data (circles), best fit QLF 
(solid line). Inset shows the fraction of the z = ma s s dens i ty (pbh(z)/pbh(0))> o n a l inear scale. Green di amond, blue square, and red cross show the local 
black hole mass density calculated in Yu & Tremaine 1 2002); Marconi et al. 1 2004), and Shankar et al. 1 2004), respectively. Bottom center shows the integrated 
black h ole mass fun ction at z = (thick black line) and z = 1, 2, 3 (blue, green, and red lines, respectively), compared to the inferred local mass function 
from Marconi et al. 12004) (yellow shaded range shows 1<t range) and Shankar et al. 12004) (thin black line). Bottom right shows the integrated cosmic X-ray 
background (black line), compared to the observations of Barcons et al. 1 2000) and Gruber et al. 1 1999) (blue and red lines, respectively, with shaded yellow lcr 
observational uncertainty). 



then the integrated black hole mass density is given by inte- 
grating the luminosity density ( Solt anll982l) . FigureflQIshows 
the result of this integration from z — > oo to z, assuming e r = 
0.1. We extrapolate beyond z = 6 but this only changes the z = 
result by ss 0.03 dex. Note that the vertical axis is logarith- 
mic; on a linear scale most of the growth oc curs at z 5, 1 ~ 2, 
in agreement with previous estimates (e.g.. | Yu & Tremaind 
2002; iMarconi et all2004t IShankar et all2004l) . 
This gives a z = black hole mass density of 

Pbh(z = 0) = 4.81!^ (— ) h 2 7l) x 10 5 M Q Mpc- 3 . (27) 

Comparing this to calculations from local bulge mass, lu- 
minosity, and v elocity dispersion functi ons: pbh = 4.2 ± 1 . 1 
JShankar et alJ l200l. q R h = 4.6+i? jMarconi et all (2004), 
Pbh = 2.9 ± 0.5 (lYu & Tremainel 120021) (each in units of 
hj x 10 5 M Q Mpc" 1 ), the agreement is good. The domi- 
nant sources of uncertainty are systematic: Tables [3] and |4] 
show the integrated Pbh(z = 0), demonstrating that differ- 
ent fitting functions with similar x 2 / l/ ( or e -g-' allowing L* 
to evolve as a polynomial in z instead of £) yields simi- 
lar factors ~ 1.5-2 systematic differences in pbhOz = 0). 
An additional contribution from Compton-thick sources is 
another uncertainty, although it cannot be much larger or 
would overpredict the present mass density and X-ray back- 
ground (bel ow). Fractional uncertainties in e r are substan- 
tial, and e.g. Hopkins, Naravan, & Hernquist (2006) estimate 
based on radiatively inefficient accretion flow models (e.g., 
lNaravan& Yl l995) a mass-weighted effective e r a factor of 
~ 0.8 lower than the "radiatively efficient" value owing to 
some growth at low Eddington ratio. Given these uncertain- 
ties, it is not necessarily meaningful to compare with the local 



black hole mass density with any greater statisti cal a ccuracy. 

Having factored out the quasar lifetime in § 14. II to deter- 
mine the rate of build-up of individual black holes 0(Mbh), 
we can trivially integrate this over redshift to determine the 
relic black hole mass function (see Hopkins et al. 2006b] for a 
full derivation of the relevant equations). We show the result- 
ing z = black hole mass function in Figure[l0| compared to 
the local mass function estimates from Shankar et al. (2004) 
and Marconi et al. (2004), determined from local bulge ve- 
locity dispersion, mass, and luminosity function observations 
(Marz keet al.irr9 94: Koc hanek et al.ll2001t ICole et alJl200U 
Naka mura et all 120031 iBlanton et all 120031 iBernardi etall 
120031 ISheth et al. 2003). The integrated black hole mass func- 
tion agrees well at all masses. We emphasize that there are 
no free parameters to fit to the black hole mass function in 
this analysis, as there are in the traditional analyses of e.g. 
lYu & T remaine (2002). The mass functions at several red- 
shifts are also shown: the highest-mass black holes are formed 
at redshifts z ~ 2, and low mass black holes assemble at z S. 1, 
in agr eement with a numb er of previous studies (e.g., [M erloni 
2004tlMarconi et al.l2004tlShankar et al.l2004tlHopkins et alJ 
2006B - 

Finally, we can integrate over the bolometric QLF, includ- 
ing our full modeling of spectral shapes and column den- 
sity distributions to determine the X-ray backgrou nd spectrum 
(see also lUeda et"aTll2003l iHoo kins et al. 2006b). Figure ITol 
compa res this with t he obs erved X-ray background spe ctrum 
from Barc ons et alJ fcOOOl) at E < lOkeV and iGruber et alJ 
(1999) at E > 3keV (normalized to the same amplitude in 
the range of overlap based on the detailed compilation in Bar- 
cons et al. 2000) with the final ~ 10% absolute normalization 
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uncertainty shown as the yellow range. The predicted X-ray 
background spectrum agrees well with that observed over the 
entire ~ 1 - 100 keV range over which "normal" AGN spec- 
tra are expected to dominate the X-ray background contribu- 
tion, and combined with the black hole mass density con- 
straint implies that the X-ray background is dominated by 
AG N with canonical radiative efficiencies e r ~ 0.1 (compare 
e.g. lElvis et alJl2002i) . This also puts strong limits on ad- 
ditional contributions of Compton-thick quasars, as a frac- 
tion much l arger than the (luminosity-dependent) ~ 30% of 
lUeda et al] ( 120031) would be problematic. 

4.3. UV and Ionizing Backgrounds 

Figure ^2 shows the specific luminosity density e„ at 
9 12 A, from our full spectral and obscuration modeling, in- 
tegrated over our full best-fit QLF model at each redshift. 
Essentially, this is identical to integrating the B-band lu- 
minosity functions and usin g the conversion L^(912A) = 



10 



18.15 



ergss-'Hz-'qWZ^) , 



(lElvis et alJ 1994; T elfer et alJ 
2002; Ric hards et alJl2006d) . We also show the bolometric 
luminosity density from Figure|S] renormalized by a constant 
"bolometric correction," and the e„ implied if this luminos- 
ity density were proportional to the number density of bright 
quasars, e„ oc $(Mg < -27), both normalized to give the in- 
ferred e„ at z = 3. The latter correction is given by 

e v « 0.95 x lO^ergss-'Hz" 1 f ^f*"? ] . (28) 

L 10~ s Mpc~ j J 

At z < 2, the faint and bright end slopes of the QLF are rela- 
tively steep, indicating an important contribution to luminos- 
ity densities from faint objects. As a result, jboi (neglecting 
obscuration) and $(Mb < -27) (neglecting all faint objects) 
are poor proxies for e u . Above z ~ 2, however, the luminos- 
ity density is dominated by bright (and relatively unobscured) 
objects, and Equation d28t provides an accurate approxima- 
tion to the full calculation of the specific luminosity density 
from all available constraints in Tabled 

Given e v at 9 1 2 A and some description of the IGM opacity, 
we can estimate the ionization rate contributed by the quasar 
population. The right panel of Figure shows this rate 



r_i2 (rate in 10 ionizations s 



the con versions from Fa rdal et alJ (Il99t 
J2001 



, where we have used 
i): lSchirber & Bullockl 



r_i 2 (z)w 2.0(1 +z) 



-1.5 



04 



3 + |auv| 



(29) 



with auv defined in § 12.11 and £24 = £;,/(10 24 ergss 'Hz '). 
We compa re with several obs e rvational estimates of the total 
F_„ fromlScott etal.l(l2000): McDonald & Miralda-Fscude 
(2001 . 120031) : iTvtler et all (120041): [Bolton et all 12005): 
iRollinde et all (120051) : iFan et all (12005^ using various meth- 
ods (for a review, see Fan et al. 2006). Note that this as- 
sumes a specific IGM optical depth, escape fraction, and Ly- 
man Limit System distribution/evolution, which are by no 
means well-constrained, so we caution that we do not mean 
this to be a rigorous calculation of T_i2 and also show the 
same estimate multiplied (divided) by a factor of 2. Ul- 
timate l^_jMs_jough£dculation_s^ previous est imates 
(e.g.. lMadau et al.ll999tlFan et al.l2007btlBolton et all2005l) 
that at z ~ 2, the ionizing background contributed by quasars 
is comparable to (~ 50%) that observed, but the quasar con- 
tribution drops much more rapidly with redshift than the ion- 
izing background. The contributed t v declines from z ~ 2 - 6 
roughly oc (1+z)" 4 (giving a factor ~ 60 drop in T-n), 
whereas the observed total T_i2 falls by only a factor ~ 10. 



5. DISCUSSION/CONCLUSIONS 

We have used a large set of observed quasar luminosity 
functions in various wavebands (Tabled, from the IR through 
optical, soft and hard X-rays, and emission line measure- 
ments, combined with recent estimates of the quasar column 
density distribution from hard X-ray and IR observations, and 
a large number of observations from the radio through hard X- 
rays determining the distribution of quasar spectral shapes to 
estimate the bolometric quasar luminosity function. This al- 
lows us to fit the bolometric luminosity function over a wide 
baseline in both luminosity and redshift, from bolometric lu- 
minosities ~ 10 41 - 10 49 ergs _1 and redshifts z = 0-6, reliably 
determining the location of the break even at z ~ 3 and con- 
straining its location with > 2er confidence up to z ~ 4.5. 

Combining observations at different wavelengths but at 
the same redshift allows us to test the self-consistency of 
measurements of the column density distribution and spec- 
tral shape. With the best-fit current estimates of the col- 
umn density distribution and quasar spectral shape, both 
of which depend on luminosity ( the unobscured fraction 
and optical/UV prominence in the quasar spectrum in- 
creasing with luminosity), a single bolometric QLF self- 
consistently reproduces the observed QLFs in each avail- 
able band at all redshifts for which we have compiled 
measurements. A constant (luminosity-independent) col- 
umn density distribution cannot self-consistently reproduce 
the observed QLF in different b ands, and is increasingly 
ruled out by direct observations (Hill. Goodrich. & DePov 
1 19961 ISimnson. Rawlings TTLacvlligggt IWillott et alJl2000t 
Simps on*& Rawlingsl 120001 ISteffen et alJ l2003t lUeda et alJ 
20031 iGrimes. Rawlings. & Wihottl 120041 iHasinge J 120041: 
Sazonov & Revnivtsev 2004; Barger et al. 2005; Simpson 
1 2005 1 : lHao et alJl2005t ISazonov et alJl2006HShinozaki et alJ 
120061 iBeckmann et afll2006bl: iBassani et alJl2006l) . A strong 
redshift dependence in the column density distribution is 
also ruled out: measu rements of such a dependence (e.g., 
La Franca et a fit the observations with compara- 

ble accuracy to a redshift-independent model up to z ~ 
1 (where the observations are generally calibrated), but 
extrapolate poorly to higher redshifts. Likewise, a sin- 
gle quasar spectrum cannot self-consistently reproduce the 
observed QLFs in different wavebands; it is necessary 
to account for the dependenc e of s p ectral shape on lu - 
minosity see n in IWilkes et alJ (11 9941): iGreen et alJ 1 1995): 
Vignal ilital] d2003h: [Strateva et alJ ll2005l) : iRichards et alJ 
(2006c); Steffe n et al.l (I2006). 

Extrapolating any fitted luminosity function outside of its 
measured luminosity and redshift range can b e inaccurate b 
orders of magnitude (see, e.g. Figure 19 of Rich ards 
(2006b)), and we have demonstrated the importance of ac- 
counting for the detailed luminosity dependence of quasar 
SEDs and obscuration. It is also important to include the scat- 
ter in bolometric corrections, i.e. the fact that there is not one 
spectral shape even at a given luminosity, as failure to account 
for this can over-predict the number density of the most lu- 
minous quasars by an order of magnitude, and over-predict 
the bolometric luminosity density by ~ 30%. Given the large 
baselines spanned by our compiled samples and our full treat- 
ment of these effects, we provide a number of fitting formu- 
lae to simplify future comparisons, including: the bolometric 
QLF itself (Table [3jl, median bolometric corrections (Equa- 
tion|2ji, the dispersion in bolometric corrections (Equation^, 
and effective "obscured/visible" fractions (Equation |3 as a 
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FIG. 1 1 . — Left: Specific luminosity density at 912 A, calculated using our full spectral and obscuration modeling (solid line). If instead we adopt a constant 
conversion from the bolometric luminosity density in Figure 151 (or. similarly, ignore obscuration) the dashed line is obtained. Taking a constant conversion 
from the number density of quasars with Mg < —27 (Equation 1281 here normalized to give the exact specific luminosity density at z = 3) yields the dotted 
line. Right: Implied r ate of i onization T_i2 (in units of 10~ 12 s -1 ) as a function of redshift, from our full modeling (black line), adopting the conversions from 
e„ from Fardal et al. 1 1998); Schirber & Bullock 1 2003). Dotted lines show t he same, but arbit rarily increased (decre ased) by a facto r of 2. Observational 
estimat es of r_i2 are shown, from Bolton et al. 1 2005) (diamonds at z = 2-41. I T vtler et al] 120041) (star at z = 1.9), Rollinde et al. 1 2005) (triangle at z = 2.75), 
McDonald & Miralda-Es cudl 1200 U I2003ir (sauares at z = 2.4 - 5.2). iFan et all 120051) (circles at z = 5 - 6), and bcott et all 1200(1) (boxes at z ~ 0- 1 and 
z ~ 2—4). Below z ~ 2, the faint quasar contribution is important, with the effects of obscuration and a changing faint end slope requiring a detailed calculation 
of the specific ionizing flux. However, above z ~ 2, the flatter bright and faint end slopes imply that bright quasars dominate the luminosity density, and the 
evolution in the number density of bright quasars becomes an accurate proxy for the evolution in e„. Quasars can account for all of the observationally estimated 
r_i2 at z < 1, and ~ 50% at z ~ 2-3 (for the specific conversions adopted), but the quasar luminosity density drops off much faster than T_i2 at higher z. 



function of luminosity. 

We find that the faint-end slope of the QLF flattens at in- 
creasing redshift, at high significance (> la). Evidence for 
such a trend comes from a number of different measurements 
JPage et alJ 119971 jMivaii et alJ 120001 120011: lLa Franca et al] 
1 20021 ICowie et alJl2003t lUeda et aljl2003t iFiore et alJl200a 
Munt et alJl2004t ICirasuolo et al jl2005t lHasinger et a!1l2005l) 
(although we still find ~ 3a evidence for this evolution with- 
out these data), but independently these samples detect such 
evolution with only marginal (< 2a) significance. The most 
signifi cant detecti o n from a single study thus far has been that 
of lHasinger et al] (|2p05), who find an improvement from a 
~ 5% 2D K-S probability to ~ 36% allowing for faint-end 
evolution; it is simply difficult to obtain the large baseline 
necessary for strong constraints from a single sample. There 
has also been some debate o ver whether the steep faint-en d 
slopes seen at z 5 , 0.2 by e.g. Sazonov & Revnivtsev (2004); 
lHao et af] i2005l) are consistent with those seen in other (par- 
ticularly optical) surveys, and indeed Figure [6] demonstrates 
they are not consistent with the PLE model. However, ac- 
counting for the evolution in the faint end slope with redshift, 
Figure [6] demonstrates that these observations are consistent 
with the others we compile. Our analysis also allows us to un- 
derstand the different shapes of the QLF in different observed 
bands, and Figures |4] [6j &0 demonstrate that with obscura- 
tion more important in low-luminosity quasars, the faint-end 
optical QLF is relatively flattened at all redshifts. Together 
with optical samples generally not probing to luminosities as 
faint as the X-ray samples, this explains why the trend of evo- 
lution in the faint-end slope has been preferentially observed 
in X-ray samples. 

This confirms, at high significance, that there is faint-end 



evolution in the data. This does not, of course, imply that the 
evolution is necessarily real, because selection effects such as 
incompleteness at low luminosities and high redshifts, or an 
evolving Compton-thick population, could play a role. How- 
ever, constraints on this quantity are interesting. For example, 
to the extent that the evolution is r eal, it supports at much 
greater significance the findings of Hopki ns et alJ (|2006a), 
namely the predicted functional dependence of the faint-end 
shape on L* (and, as a consequence, on redshift). This model 
is motivated by hydrodynamic simulations (Hopkins et al. 
2005a,b,c) and analytical models of accretion rate evolution in 
feedback-driven outflows (Hopkins & Hernquist 2006), which 
predict both that the faint-end shape should be controlled by 
the shape of the quasar light curve, itself a function of quasar 
and host galaxy properties. To the extent that the faint-end 
slope does or does not evolve, it directly constrains the quasar 
"lifetime" as a function of luminosity (i.e. differential time 
quas ars spend per l uminosity interval in both growth and de- 
cay; [Hopkins et al. 2005c), and consequently feedback cou- 
pling modes and their impact on the host galaxy. 

We also find at high significance (> 6a) that the bright end 
slope of the QLF evolves with redshift, becoming shallower 
at high redshifts z > 3. This agrees with the early suggestions 
of lSchmidt et alJ dl995l) : lFan et alJfcOOlalbl). and s ubsequent 
more robust measurements of lpnrhardTetal] fc006bl) . We in- 
fer significant (> 3a) evidence for such evolution even with 
those samples removed (& 7ll i = 1.22±0.15, k lu 2 =— 0.80± 
0.09), but in this case the result is heavily influenced by the 
X-ray data, which at the highest luminosities has a large bolo- 
metric correction (L\, \ ~ IQQLhx, from the observed a ox -L 
relation). Since the QLF is falling steeply at these luminosi- 
ties, a small change in the bolometric correction will signifi- 
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cantly change the inferred slope; Figure |4] aptly demonstrates 
that the optical and IR samples provide a much more robust 
probe of the "true" bright end slope. Note that the best-fit 
trend at z ~ 6 favors a sha llower slop e than recently sug- 
gested bv lFan et al] (2004): C ool et ali (120061) . although the 
small number of objects make this only a ~ la effect. Further 
wide-area surveys are needed to confirm the trend. 

There is also marginal evidence that the bright-end QLF 
slope is somewhat shallow at z ~ 0, steepening to z ~ 0.5 - 1, 
then becoming shallower with z as describe d above. This 
low-red shift trend is marginally evident in the Richar ds et alJ 
(2006b) sample as well (see Figure^, although with a smaller 
low-z baseline the authors assumed a constant bright-end 
slope up to z ~ 3. There are a number of potential system- 
atic effects which complicate understanding this trend. For 
example, imperfect accounting of the (potentially quite rapid) 
QLF evolution over a narrow redshift interval will bias the lu- 
minosity function to a higher number density at large L, as 
will many general binning procedures used in presenting the 
data, as well as host galaxy contamination (see e.g. Figure 1 1 
of Hao et al. 2005). These problems are well known from 
local galaxy luminosity functions, and can introduce signifi- 
cant effe cts even in samples s panning the range e.g. z ~ 0-0. 1 
(see e.g.. lBlanton et alJ200 3). where many of the low-redshift 
QLF measurements cover a significantly larger z ~ 0-0.3. 
Furthermo re, the fain t-end slope is observed to be steep at 
z ~ (e.g. JHao et al. 2005), implying a weak break and large 
covariance between the best-fit parameters (much larger than 
at any other redshifts where the break is constrained). More 
careful re-analysis of t hese samples and e.g. recent local lumi- 
nosity func tions from | Sh inozaki et alJ ( 12006ft . ISazonov et alJ 
(2006), and Beckmann et al. (2006b), accounting for these ef- 
fects, is needed. 

Again, despite systematic uncertainties, evolution in the 
bright-end shape is interesting. Models predict that it should 
evolve wit h redshift, as a consequenc e of feedback from AGN 
(e.g.. iKauffmann & Haehiieltl 120001 IWvithe & Loebl 120031: 
IScannapieco et alJl2005l) . Athigh redshift, with small galax- 
ies, large gas content, and rapid merger rates, the quasar ac- 
tivity is expected to trace dark matter halo growth and merger 
rates, until a redshift z ~ 2 - 3 when feedback from quasar 
growth or other mechanisms begins to dominate. This "shuts 
down" activity, termin ates subsequent black h ole growth and 
local s tar formation ( Hopkins e^ alJ |2QQ6d |$prin gel et al~l 
12005 aft , injects entropy (Scannapie co & OH I2004) and radio- 
mode feedback fro m low accretion rate "dying" black holes 
llCroton et alJl20 06) heats the surrounding IGM and prevents 
new infall, suppressing further galaxy and BH growth. In 
these models, this occurs first in the largest systems, which 
evolve the most rapidly and have the most massive, vio- 
lent black holes, leaving a gas-rich merger and correspond- 
ing qu asar history which shi fts to lower masses at lower red- 
shifts (Hopkins et al. 2006d), and building up lower-mass el- 
liptical ga laxies at preferentially later times ( Hopki ns et alJ 
2006c aa), explaining observ ed trends of "cosmic downsiz- 
ing" ("e.g.. ICowie et alJl!996ft in the context of hierarchical 
structure formation. 

Our estimate of the bolometric QLF allows us to imme- 
diately derive several interesting quantities. The bolometric 
quasar luminosity density (and correspondingly, total black 
hole accretion rate) of the Universe follows a similar qualita- 
tive trend to that of L„, peaking at a well-defined z ~ 2.154 ± 
0.052. The integrated black hole mass density to z = is 
consistent with local estimates from bulge mass and velocity 



dispe rsion functions (e.g.. lMarconi et all2004tlShankar et alJ 
2004), with room for observationally estimated Compton- 
thick fractions of ~ 30-40% which we have included (but 
not, e.g. > 60%). The predicted X-ray background spectrum 
agrees well with t hat observed, in a similar manner to that in 
lUeda et all ll2003ft . The UV luminosity density and ioniza- 
tion rate are consistent with all existin g constraints, and ac- 
count for observed ionization rates from lScott et alJ (|2()00) at 
Z < 1 and a significant (~ 50%) fraction at z ~ 2-3, but subse- 
quently dropping more rapidly than observed rates. Because 
of the flattening of the bright and faint-end QLF slopes with 
redshift, the luminosity density at high redshift is dominated 
by bright quasars, and above z ~ 2.5 the number density of 
bright, optical quasars (e.g. Mb < -27) becomes an accurate 
proxy for the UV lumino sity density, related by a conversion 
which we calibrate (see § I4.3> . 

We can also de-convolve the bolometric QLF with the 
quasar lifetime as a function of lumin osity and black hole 
mass; similar to the procedure in e.g. lYu & Tremainel (2002) 
but adopting the physic ally motivated quasar lifetimes from 
Hopkins et al. (2006a b) to determine the "birthrate" (forma- 
tion rate) of black holes as a function of mass and redshift. 
The resulting evolution of the black hole mass function illus- 
trates the feedback-driven scenario above, with rapid buildup 
to z ~ 3, then activity shutting down in the highest mass sys- 
tems, and the characteristic mass of black holes "in forma- 
tion" shifting from ~ 1O 9 M at z ~ 2 to — 2 x 10 7 M Q at 
z ~ 0. The resulting bla ck hole mass function agrees well 
with t hat inferred locally (Ma rconi et al]l20"04llShankar et alJ 
2004). 

Measurements of the QLF are converging towards a self- 
consistent picture of the QLF shape and evolution across ~ 8 
orders of magnitude in luminosity and ~ 9 in space density, 
frequencies from ~ 24 /im to ~ lOOkeV, through mid and 
near-IR, optical, UV, soft and hard X-rays. Except at the 
highest redshifts z > 4, where statistics are still poor and faint 
luminosities have not been probed, systematic errors will al- 
most certainly overwhelm statistical errors. Indeed, compar- 
ison of samples with overlapping redshift and luminosity in- 
tervals shows that while the agreement is good, the systematic 
offsets (albeit small) between the largest samples are much 
larger than their typical statistical errors. Some of these ef- 
fects, such as cosmic variance, are at least straightforward 
to understand and reduce. However, others are more funda- 
mental. It may not be possible to define a bolometric QLF to 
any greater accuracy than we have done (again, excepting the 
highest redshifts where larger samples would improve things 
considerably), at least in the near future. 

For example, on both small scales characteristic of accre- 
tion disks, black hole driven winds, and molecular outflows 
or tori, and large scales characteristic of potential merger- 
driving and galaxy-scale quasar obscuration, quasars are not 
isotropic. So even a complete spectrum does not imply the 
bolometric luminosity, lacking a more complete model for the 
spectrum as a function of viewing angle and other quantities 
(e.g. accretion rate, mass, spin). Perhaps by considering the 
distribution of luminosities and accounting for the fact that 
there is not one spectral shape for a given total bolometric lu- 
minosity, we have statistically accounted for some of this, but 
it is by no means clear, and will not necessarily be so if there 
are correlations between e.g. the column density distribution 
and intrinsic Eddington ratio or spectral shape distribution, or 
if e.g. the dependence of spectral shape on luminosity is re- 
ally driven by an accretion rate or viewing angle dependence. 
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These are important but difficult and often degenerate ques- 
tions, which can only be answered with more detailed obser- 
vations of quasar spectra at many different frequencies, cou- 
pled with detailed modeling of accretion processes and quasar 
fueling. 

Finally, in order to enable simple comparison of our results 
with future theoretical models and observations at a wide va- 
riety of wavelengths and redshifts, we provide for public use 4 
a simple "QLF calculator" script to return the QLF at a given 
redshift, in an arbitrary observed band or frequency. This cal- 
culates the observed QLF using the full modeling of bolo- 
metric corrections and extinction discussed herein, for any of 
the QLF fits in Tables [3] & E] or for an arbitrary bolometric 
QLF. As we have emphasized for the large number of histori- 
cal fits to the QLF, the fits therein should not be extrapolated 

http : / /www . cf a . harvard . edu/ -phopkins /Site/qlf . html 



beyond the observed ranged in Figures®&\7\ for example to 
very high redshifts (z > 6), low luminosities at high redshifts 
(z > 4.5), and unconstrained wavelength regimes such as the 
extreme IR, UV, or hard X-rays. 
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TABLE 1 
Measurements of the QLF 



Reference 


Survey/Field" 


Rest Wavelength/Band 


z Range b 


Luminosity Range b 




ZVagn 


Plotting Symbol 


Optical: 


Cristiani et al. (2004) 


GOODS 


1450 A 


~4-5.2 


-21 > M 1450 > -23.5 


0.58/1 


1-4 


crosses 


Croom et al. (2004) 


2QZ/6QZ 


B 


0.4-2.1 


-20.5 > M g > -28.5 


23.1/10 


20,905 


asterisks 


Fan et al. (2001a) 


SDSS (Equatorial Stripe) 


1450 A 


3.6-5.0 


-25.5 > M 145 o > -27.5 


6.21/9 


39 


pentagons 


Fan et al. (2001b, 2003, 2004) 


SDSS (Main & Southern Survey) 


1450 A 


~ 5.7-6.4 


-26.5 > M 145 o > -28 


2.12/3 


9 




Hunt et al. (2004) 


LBG survey 


1450 A 


~2-4 


-21 > M 145 o > -27 


4.74/6 


11 


diamonds 


Kennefick et al. (1995) 


POSS 


B 


4.0-4.5 


-26.5 > M B > -28.5 


14.8/2 


10 


triangles 


Richards et al. (2005) 


2dF-SDSS 


g 


0.3-2.2 


-21 > M g > -27 


137/99 


5,645 


circles 


Richards et al. (2006b) 


SDSS (DR3) 


i'(z = 2)~2500A 


0.3-5.0 


-22.5 > Mi > -29 


247/101 


15,343 


squares 


Schmidt et al. (1995) 


PTGS 


B 


~ 3.5-4.5 


-25.5 >M B > -27.5 


8.04/4 


8 


inverted triangles 


Siana et al. (2006) 


SWIRE (ELIAS-N1/N2) 


1450 A 


~ 2.8-3.4 


-23.5 > M 1450 > -26.5 


4.74/6 


~100 


crosses 


Wolfetal. (2003) 


COMBO- 17 


1450 A 


1.2-4.8 


-23.5 > M 1450 > -28.5 


54.2/27 


192 


stars 



Soft X-ray: 



Hasinger et al. (2005) 


ROSAT (RASS+RDS) + CDF-N/S 


0.5 


-2keV 


0.015 


-4.8 


10 42 <Lo.5- 


.2 < 10 48 ergs-' 


169/51 


2,566 


circles 


Miyajietal. (2000,2001) 


ROSAT (RASS+RDS) 


0.5 


-2keV 


0.015 


-4.8 


10 41 < Lo. 5 _ 


, 2 < lO^ergs" 1 


112/41 


691 


stars 


Silverman et al. (2005b) 


CHAMP+ROSAT (RASS) 


0.5 


-2keV 


0.1- 


-5 


10 44 - 5 < Lo.s 


_ 2 < 10 46 ergs-' 


24.1/9 


217 


squares 



Hard X-ray: 



Barger et al. (2003a,b) 


CDF-N 


2- 


-8keV 


~5- 


-6.5 


10 43 <L 2 -8 < lO^ergs" 1 


1.02/1 


1 


diamonds 


Barger et al. (2005) 


CDF-N/S + CLASXS + ASCA 


2- 


-8keV 


~0.1 


-1.2 


10 42 <L 2 _ 8 < lO^ergs- 1 


41.0/30 


601 


squares 




CDF-N/S + CLASXS 


2- 


-8keV 


~ 1.5 


-5.0 


10 42 <L 2 -8 < lO^ergs- 1 


15.5/9 


~100 




Barger & Cowie (2005) 


CDF-N/GOODS-N 


2- 


-8keV 


~2 


-3 


10 43 <L 2 -8 < 10 44 - 5 ergs-' 


1.73/1 


136 




La Franca et al. (2005) 


HELLAS2XMM 


2- 


lOkeV 


0.0- 


4.0 


10 42 <L 2 -io < 10 465 ergs- 1 


14.4/18 


508 


stars 


Nandra et al. (2005) 


GWS + HDF-N 


2- 


lOkeV 


2.7- 


3.2 


10 43 <L 2 -io < 10 44 - 5 ergs-' 


0.77/1 


15 


crosses 


Sazonov & Revnivtsev (2004) 


RXTE 


3- 


20keV 


0.0- 


0.1 


10 41 <L 3 - 20 < 10 46 ergs-' 


9.75/10 


77 


inverted triangles 


Silverman et al. (20()5a,c) 


CHAMP 


0.3- 


-8.0keV 


0.2- 


4.0 


10 42 <L .3-8 < 10 45 - 5 ergs' 1 


26.3/15 


368 


triangles 


Ueda et al. (2003) 


HEAO 1 + AMSS-n/s + ALSS 


2- 


lOkeV 


0.015 


-3.0 


10 41 ' 5 <L 2 -io < lO^ 'ergs- 1 


26.5/35 


247 


circles 



+ ASCA + CDF-N 



Mid-IR: 



Brown et al. (2006) Spitzer Bootes (NDWFS) 8 /an ~l-5 10 45 < L 8/am < 10 47 ergs" 1 3.77/10 183 circles 

Matute et al. (2006) RMS + ELIAS + HDF-N/S 15/mi ~ 0.1-1.2 10 42 < L 15Mm < 10 47 ergs" 1 23.4/18 148 squares 

Emission Lines: 



Hao et al. (2005) SDSS (main galaxy sample) Ha - 0.33 1 5 < L Ha < 10 9 L Q 29.5/21 ~3000 pentagons 

[Oil] ... 10 5 <Lon < W S L Q 

[OIII] ... 10 5 <Lora < 10 9 L© 



"For a detailed description of each sample, we direct the reader to the listed references (and references therein). 

b Redshift and luminosity ranges listed are for the entire sample in each case, they should not be taken to imply that the observations simultaneously span both ranges. 
c Reduced \ 2 of binned QLF with respect to our full best-fit. 



TABLE 2 

Best-Fit QLF at Various Redshifts 



w 






7l 


72 


->/ 2 I,, 
X jv 


0.1 


-5 .45 ±0.28 


11. 94 ±0.21 


0.868 ±0.050 


1.97±0.17 


89/73 


0.5 


-4.66 ±0.26 


12.24±0.18 


0.600 ±0.136 


2.26 ±0.23 


124/66 


1.0 


-4.63 ±0.15 


12.59±0.11 


0.412±0.122 


2.23 ±0.15 


182/69 


1.5 


-4.75 ±0.19 


12.89±0.13 


0.443 ±0.145 


2.29 ± 0.20 


214/86 


2.0 


^.83 ±0.05 


13.10±0.04 


0.320 ±0.046 


2.39 ±0.07 


66/67 


2.5 


-4.96 ±0.14 


13.13±0.09 


0.302 ±0.091 


2.30±0.15 


72/53 


3.0 


-5.23 ±0.12 


13.17±0.10 


0.395 ±0.060 


2.10±0.12 


45/53 


4.0 


-4.66 ±0.37 


12.39±0.32 


-0.254 ±0.736 


1.69±0.18 


54/32 


5.0 


-5.38 ±1.19 


12.46± 1.10 


0.497 ±0.458 


1.57 ±0.41 


14/13 


6.0 


-5.13±0.38 


11.0 


0.0 


1.11 ±0.13 


5/3 



a Mpc 3 

b L Q =3.9x lO^ergs" 1 



TABLE 3 
Best-Fit QLF to all Redshifts 



Model 


log</>* a 


(logL*)o b 


k L A 


kL.2 


h,3 


(71 )o 




(72)0 


£72,1 


£ 72 , 2 


Pbh(z = 0) c 


X 2 /" 


PLE 


-4.733 ±0.101 


12.965 ±0.074 


0.749 ±0.084 


-8.03 ±0.35 


-4.40 ± 1.05 


0.517 ±0.065 





2.096 ±0.083 








5.66 


1924/511 


Faint 


-4.630 ± 0.075 


12.892 ±0.058 


0.717 ±0.069 


-8.10±0.29 


-3.90 ±0.85 


0.272 ±0.073 


-0.972 ±0.268 


2.048 ±0.063 








4.97 


1422/510 


Bright 


-4.930 ±0.070 


13.131 ±0.051 


0.360 ±0.095 


-11.63 ±0.48 


-10.68 ± 1.04 


0.605 ± 0.044 





2.350 ±0.087 


1.53 ±0.14 


-0.745 ±0.081 


5.47 


1312/509 


Scatter d 


-4.815 ±0.040 


13.064 ±0.030 


0.356 ±0.061 


-11.69 ±0.27 


-9.18±0.58 


0.351 ±0.032 


-0.826 ±0.092 


2.359 ±0.052 


1.534 ±0.070 


-0.889 ±0.042 


4.63 


445/508 


Schechter 13 


-3.579 ±0.092 


11.482±0.162 


-1.78 ±0.29 


-21.22 ±1.62 


-25.93 ±3.08 


0.013 ±0.011 


-3.760 ±0.859 


0.354 ±0.020 


1.794 ±0.087 


-0.784 ±0.063 


1.16 


1254/507 


Full 


-4.825 ±0.060 


13.036 ±0.043 


0.632 ±0.077 


-11.76±0.38 


-14.25 ±0.80 


0.417 ±0.055 


-0.623 ±0.132 


2. 174 ±0.055 


1.460 ±0.096 


-0.793 ±0.057 


4.81 


1007/508 



a Mp<T 3 

b L Q = 3.9 X 10 33 ergs-' 
c lO 5 M MpiT 3 

d This fit allows for an intrinsic sample-to-sample systematic normalization variance of 0.05 dex. However, this was applied uniformly to the samples of Table HI and under-weights much of the most well-constrained 
data in comparison to the "Full" fit, which has a higher y 2 / u but more faithfully represents those data sets. 
e Adopting the modified Schechter function fit of Equation 1211 instead of a double power law. 

The resulting observed QLF at arbitrary frequency and redshift from any of these fits can be calculated using the QLF calculator script available for download at 

http : / /www . cf a . harvard . edu/~phopkins /Site/ qlf . html 



TABLE 4 
Best-Fit LDDE QLF 



Model 


log0* a 


logL„ b 


71 


72 


logL t b 


ZcO 


a 


PU6 


p2 46 


ft 


Pi 


Pbh(z = 0) c 


X> 


PDE 
LDDE 


-6.66 ±0.28 
-6.20 ±0.15 


46.64 ±0.19 
45.99±0.10 


0.858 ±0.075 
0.933 ±0.045 


2.09 ± 0.20 
2.20±0.14 


46 

46.72 ±0.05 


1.852 ±0.068 
1.852 ±0.025 




0.274 ±0.025 


4.13±0.29 
5.95 ±0.23 


-2.53 ±0.46 
-1.65 ±0.21 




0.29 ±0.34 




-0.62 ±0.17 


9.78 
5.09 


3255/511 
1389/507 



a Mpc- 3 
b ergs~' 

c 10 5 Mq Mpc~ 3 . In this formulation, the luminosity density is nearly divergent at low L, so the integration is truncated at 10 41 ergs 



The resulting observed QLF at arbitrary frequency and redshift from any of these fits can be calculated using the QLF calculator script available for download at 

http : / /www . cf a . harvard . edu/~phopkins/Site/qlf . html 



TABLE 5 
Best-Fit Quasar Formation Rate 



Model 


(log<Mo a 




(logM„) b 




^M,2 


m c 






k Vl .2 


PBHfe = 0) d 


X 2 /" 


Fixed 72 
Variable 72 
Schechter 6 


-3.808 ±0.037 
-3.830 ±0.031 
-3.571 ±0.080 


-4.46 ±0.31 
-4.02 ±0.36 
-3.81 ±0.46 


8.894 ±0.037 
8.959 ±0.032 
8.710±0.130 


1.54±0.12 
1.18 ±0.13 
1.33±0.21 


-4.46 ±0.31 
-6.68 ± 0.44 
-7.64 ±0.99 


0.20 
0.20 
0.20 


2.44 ±0.11 

2.86±0.16 
0.79 ±0.10 




1.80±0.18 
1.33±0.19 




-1.13 ±0.09 
-1.17±0.12 


6.19 
6.28 
4.27 


1338/511 
1206/509 
1434/509 



"Mpc" 3 Gyr" 1 
b M Q 

c In this model, the QLF faint-end slope is dominated by the shape of the quasar lightcurve, giving only weak constraints on rji . We adopt the maximum acceptable rj\ 0.20 
(lcr) as an upper limit, and show the fit results for this case. 
d 10 5 Mo Mpc~ 3 

c Adopting the modified Schechter function fit of Equation l2 1 l instead of a double power law. 



